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ABSTRACT 

We present analytic calculations of angular momentum transport and gas inflow in galax- 
ies, from scales of ~ kpc to deep inside the potential of a central massive black hole (BH). 
We compare these analytic calculations to numerical simulations and use them to develop a 
sub-grid model of BH growth that can be incorporated into semi-analytic calculations or cos- 
mological simulations. Motivated by both analytic calculations and simulations of gas inflow 
in galactic nuclei, we argue that the strongest torque on gas arises when non-axisymmetric 
perturbations to the stellar gravitational potential produces orbit crossings and shocks in the 
gas. This is true both at large radii ~ 0.01 — 1 kpc, where bar-like stellar modes dominate the 
non-axisymmetric potential, and at smaller radii < 10 pc, where a lopsided/eccentric stellar 
disk dominates. The traditional orbit crossing criterion is not always adequate to predict the 
locations of, and inflow due to, shocks in gas+stellar disks with finite sound speeds. We de- 
rive a modified criterion that predicts the presence of shocks in stellar dominated systems even 
absent formal orbit crossing. We then derive analytic expressions for the loss of angular mo- 
mentum and the resulting gas inflow rates in the presence of such shocks. We test our analytic 
predictions using hydrodynamic simulations at a range of galactic scales, and show that they 
successfully predict the mass inflow rates and quasi-steady gas surface densities with a small 
scatter ~ 0.3 dex. We use our analytic results to construct a new estimate of the BH accretion 
rate given galaxy properties at larger radii, for use in galaxy and cosmological simulations and 
semi-analytic models. While highly simplified, this accretion rate predictor captures the key 
scalings in the numerical simulations. By contrast, alternate estimates such as the local vis- 
cous accretion rate or the spherical Bondi rate fail systematically to reproduce the simulations 
and have significantly larger scatter. 
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. ^ 1 INTRODUCTION 

^N Gas inflows into galactic nuclei play a critical role in galaxy forma- 
jrt tion, initiating a wide range of phenomena including starbursts, the 
formation of bulges, spheroidal galaxies, and nuclear star clusters, 
and accretion onto supermassive black holes (BHs). The initial trig- 
gers for these inflows are diverse and may include galaxy mergers 
and tidal encounters; bar and spiral disk instabilities; secondary and 
tertiary instabilities ("bars within bars"); and large-scale triaxial or 
binary structures in galaxy halos. 

Intense star formation in galactic nuclei also has a major effect 
on the structure of the resulting systems, with dissipative processes 
dominating the formation of the inner ~kpc of bulges at surface 



mass densities reaching ~ 10 Mq kpc (Ostriker 



1980 



Carl- 



[bergl[T986l |Gunnl[T987l |Kormendy|[T989"l |Hernquist et al.||J993 
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Hopki ns" et al.||2009a|c) . There is also strong indirect evidence 
that the most luminous starbursts, ultra-luminous infrared galaxies 
(ULIRGs), are accompanied or followed by periods of luminous 
quasar activity (San ders et al.|19 88a b, Dasyra et al. 2007). 

On galactic scales, numerical simulations have shown that in 
gas-rich mergers, tidal torques lead to rapid gas inflow into the 
central < kpc l |Hernquist|1989"l|Barnes & Hernquist|1991|[T996l >. 
Minor mergers and instabilities in self-gravitating disks (bars and 
spiral waves) can drive similar inflows (Hernquist & Mihos 1995; 
Bourn aud et al.||2005[ |Eliche-Moral et al. 2008, Younge r et al.| 
2008). And in high redshift, gas-rich disk galaxies, the sinking of 
massive star-forming clumps may produce analogous strong gravi- 
tational torques (Bournaud et al. 2007, Genzel et al. 2008). 

How gas is transported from ~ 1 kpc to the much smaller 
scales relevant for massive BH growth remains uncertain (e.g., 
Goodman 2003 ). Once gas reaches sub-kpc scales, the torques from 
large-scale bars, and even from major mergers, become less ef- 
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ficient. Local viscous stresses - which are believed to dominate 
angular momentum transport near the central BH (e.g., Balbus & 



Hawley 19 98} - are ineffi cient at radii > .01 — 0.1 pc (e.g.,|Shlos 
man & Begelman 1989) |Goodman||2003| [Thompson et al.|2005) . 



And scattering of individual gas clumps or molecular clouds cannot 
produce the high accretion rates needed to explain bright quasars 
(Hopkins & Hernquist|2006[|Kawakatu & Wada|2008||Nayakshin 
& King 2007 1. As a consequence, most models implicitly or explic- 
itly invoke some form of global gravitational torques ("bars within 
bars"; Shlosman et al. 19891 to continue transport to smaller radii. 
Recently, numerical simulations using a variety of techniques have 
been able to follow the inflow of gas to smaller scales in galaxies 



Ipscala et al.||2004| |Escala|[2007{ |Mayer et al.|[2007l |Wise et al.| 



|2008[|Levine et al.|2008||Hopkins & Quataert|2010a| [ This has al- 
lowed simulations to begin to address some of the long-standing 
questions about gas inflow in galactic nuclei. 



In a previous paper (Hopkins & Quataert 2010a), we presented 
a large suite of ~ 100 simulations that begin from initial condi- 
tions of either gas-rich major mergers or bar-unstable disks, and 
follow the inflow of gas from > lOkpc to < 0.1 pc; at the small- 
est radii, the system begins to resemble a traditional, viscous ac- 
cretion disk. Our method uses a series of independent numerical 
simulations on successively smaller scales (rather than an AMR 
or particle-splitting approach), which allows us to systematically 
survey how a variety of initial conditions and galaxy properties af- 
fect gas dynamics in galactic nuclei. These calculations show that 
large-scale bars can indeed drive gas to <kpc; and if the disky ma- 
terial at smaller radii becomes self-gravitating, it triggers a cascade 
of secondary gravitational instabilities that transports gas deeper 
in the potential well. This qualitatively resembles the "bars within 
bars" scenario. However, there are several critical differences. First, 
the secondary instabilities on intermediate scales (~ 10— 100 pc) 
exhibit a diverse range of morphologies. A more accurate charac- 
terization would thus be "it's non-axisymmetric features all the 
way down" (or "stuff within stuff"). This is consistent with the 
fact that observed systems that do show nuclear asymmetries can 
show nuclear rings, barred rings, one or three-armed modes, and 
clumpy /irregular structures (Martini & Pogge 1999. Peletier et al. 



1999} [K napen et al.|2000| [Laine et al.|2002| [Knapen et al.|2002| 
Green e et al.|2009f . Moreover, in our simulations, the inflow rates 
and gas morphology are time-variable, leading to modest duty cy- 
cles that may explain the generally weak correlation in observa- 
tions between mode structure and inflow/outflow on different spa- 
tial scales. A second key difference between our numerical results 
and the canonical "bars within bars" model is that once gas reaches 
the radii where the gravity of the BH begins to dominate (~ 10 pc), 
bar-like modes no longer survive. However, a new instability ap- 
pears, in the form of a slowly precessing lopsided/eccentric disk or 
one-armed spiral (m = 1) mode, which dominates the angular mo- 
mentum transport to < 0.1 pc. The resulting torques can drive gas 
inflows at rates comparable to those required to sustain the brightest 
observed quasars, ~ IOMq yr _1 . And the surviving stellar relics 
of these modes closely resemble the observed nuclear disks in M3 1 
(Lauer et al. 1993 1 and other galaxies (Hopkins & Quataert 2010b). 

Given the simplifying assumptions necessary in numerical 
simulations (e.g., sub-resolution star formation and feedback mod- 
els), it is clearly important to understand analytically the physics of 
angular momentum transport. Moreover it is not possible to sim- 
ulate the entire parameter space relevant to galactic inflows, so a 
way of generalizing the numerical results would be useful, both for 
gaining physical intuition and for developing subgrid models that 



can be used in calculations that lack the physics/resolution to di- 
rectly capture gas inflow on sub-kpc scales. 

There is at present no analytic theory that accurately predicts 
the inflow rates of gas in mixed stellar-plus-gas systems, despite the 
rich history of analytic studies of gravitational instabilities in disks 
(see e.g. Goldreich & Lynden-Bell 1965a b ; Lin et al. 1969 ; Toomre 



)re 
50; 



1969; Lynden-Bell & Kalnajs 1972, Goldreich & Tremaine 1980 
Binney & Tremaine|1987||Shu et al.|1990||Toomre|1977||1981[|Lau 1 
& Bertin|1978||Adams et al.|1989[|shu et al.|1990||Ostriker et ~ 



1992). Almost all previous calculations have simplified the problem 



dramatically by focusing on strictly stellar or gas disks (as opposed 
to gas+stellar disks) and by ignoring the presence of shocks and 
cooling/dissipation in the gas. In addition, the majority of the an- 
alytic development and intuition has focused on the case of weak 
perturbations. In these limits, the leading-order instantaneous grav- 
itational torques associated with the asymmetric gravitational po- 
tential cancel exactly. Gravitational angular momentum exchange 
occurs in narrow resonances (Goldreich & Tremaine 19791, and 
the net rate at which the gas loses angular momentum is second- 
order in the fractional magnitude of the perturbation (denoted \a\ 
throughout; Lynden-Bell & Kalnajs 1972). It is typically further 
suppressed by strong powers in the radial and azimuthal wavenum- 
bers of the perturbation (in both the global and local perturbation 
limits; Kalnajs 1971). These results lead to the prediction of rel- 
atively weak torques, with sub-sonic mean gas inflow velocities 
Vinfiow < 0. 1 \a\ 2 <2~' c, (where Q is the standard Toomre Q param- 
eter and c s <C V c is the gas sound speed; Goodman 2003). 

In the numerical simulations described above of galaxy merg- 
ers, disk instabilities, and secondary and tertiary instabilities, none 
of the simplifications often utilized in the study of angular mo- 
mentum transport are in fact realized. Most of the angular momen- 
tum exchange occurs in the non-linear regime in which the non- 
axisymmetric perturbations to the gravitational potential due to the 
stars are sufficiently large that they induce orbit crossing, shocks, 
and dissipation in the gas. The multi-component dissipative nature 
of the system is thus critical. In, for example, instabilities of galac- 
tic disks, the gas orbits respond to the asymmetric potential of the 
collisionless (stars plus dark matter) material; when the asymme- 
try is sufficiently large, the gas is driven into shocks, which dissi- 
pates energy and allows the torques from the asymmetric potential 
to remove angular momentum efficiently. The "pile up" of mate- 
rial in these shocks manifests itself observationally as dust lanes 
and enhanced star formation in the leading edge of spiral arms 
l Roberts et al.|1979[[Alhanassoula|1992bl|Wada|l"994l|Bournaudl 



|et al.|20 05 ). The same basic process also dominates the angular mo 
mentum loss in major galaxy mergers, with the asymmetry in the 
stellar potential being excited by the companion (Hernquist 1989; 
|Barnes & Hemquist|199l1[l996")|Hopkins et al.|2009b) . The result- 
ing mean inflow velocities are typically super-sonic, with inflow 
rates up to ~ 100 — 1000 times larger than the predictions based 
on linear torques in single-component systems! The simulations 
also explicitly demonstrate that "hydrodynamic" pressure forces do 
not contribute significantly to the torques and inflow (see Barnes & 
Hernquist|1996||Berentzen et al.|2007J|Hopkins et al.|2009b[|Hop 



kins & Quataert 2010a}: the process driving angular momentum 
loss is overwhelmingly gravitational. 

The basic processes summarized above - gas responding 
to gravitational torques predominantly from stars, "piling up" in 
shocks, and exhibiting supersonic mean radial inflow/outflow ve- 
locities as a result - have been directly observed in a variety of 
nearby galaxies, where the role of gravitational torques can be 
mapped ( |Haan et al.|2009{|Durbala et al.|2009[|Block et al.|2001[ 
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Sakamoto et al. 1999, Quill en et al.|1995) . These observations also 
extend to galactic nuclei where secondary instabilities ("bars within 
bars") may become important ( Laine et al. 2002 , Schinner er et al.| 
|2006| Knapen et al.|2002 1, and to the gas surrounding bright AGN 
(Combes e t al.|2004| Garcia-Bu rillo et al.|2 005). The torques in- 
ferred in many of these systems correspond to the gas losing a 
large fraction of its angular momentum in a few dynamical times. 
This is much stronger than can be accounted for by torques in lin- 
ear, single-component systems and highlights the need for a better 
physical understanding of gravitational torques and inflow under 
realistic galactic conditions. 

In this paper, we develop an analytic model that attempts to 
account for the realistic complications of angular momentum re- 
distribution in systems with both collisional (gas) and collisionless 
(stellar/dark matter) components, in the presence of shocks and dis- 
sipation. This is an analytic companion to the numerical simula- 
tions in Hopkins & Quataert (2010a), which in turn builds on the 
large history of analytic and numerical studies of angular momen- 
tum transport in galaxies described above. The remainder of this 
paper is organized as follows. In §[2] we calculate the response of 
gas to general asymmetric perturbations, deriving the conditions for 
shocks and the resulting angular momentum transport; we consider 
the linear regime of no shocks (§ \--\ . strong orbit crossings that 
induce shocks (§ |2.3fr , and marginal orbit crossings (§ |2.3.3| l. We 
then briefly describe the balance between inflow and star formation 
in gaseous disks, and the resulting quasi-steady gas density distri- 
bution (§ 13). In § H] we describe our suite of high-resolution hy- 
drodynamic simulations of mergers and unstable disks, with multi- 
scale re-simulations of secondary instabilities, which we use to test 
our analytic scalings as a function of galaxy properties and phys- 
ical scale. In §[5] we use the results of §[2] to develop an analytic 
model that predicts the inflow rate to the BH given the conditions 
at large radii in a galaxy; we show that this predictor is much more 
accurate than the spherical accretion prediction or a local viscous 
model. Finally, we summarize and discuss our results in §[6] 

2 GAS RESPONSE IN ASYMMETRIC SYSTEMS 

2.1 Definitions 

We adopt a cylindrical coordinate frame, (R, 0, z), and con- 
sider an initially axisymmetric disk with an arbitrary spherical 
(BH+bulge+halo) component. We use the test particle approxi- 
mation to calculate the response of the disk to a specified non- 
axisymmetric perturbation to the potential, largely neglecting pres- 
sure forces. Since we are interested in the behavior of cold, ro- 
tationally supported, gas, the gas of interest is always in a thin 
(h/R < 1) disk, aligned with the stellar disk. 

The initial potential in the disk plane can be written $o = 
$o(R), and other properties are defined using standard notation: 



ft = 



<9$ 



R 



GM mc (<R) 



R 



'dyn 



= V C /R 

2 



■.£ + «*-£+» rf 



(i) 

(2) 
(3) 
(4) 



where V c is the circular velocity, ft the angular velocity, and k the 
epicyclic frequency. We now consider a global perturbation in the 
disk plane and define the perturbed potential 



In what follows, the subscript "1" will denote all perturbed quanti- 
ties. We will consider a frame rotating with the perturbation pattern 
speed ftp, so that in these coordinates the perturbation is stationary. 
An unperturbed orbit in the rotating frame is at fixed R with angular 
speed 0o = ft — ftp. We can represent the perturbed potential as a 
sum over modes: 



$i = a(R) $o/(0) = $o(R) J2 a ">( R ) cos ( m <; 



(6) 



For all perturbed quantities X, we define X\(R,<f>) — 
X a (R) cos (m(j>). Thus <!>,, = |a| $o, etc. Note that for a given 
particle, the perturbed and unperturbed quantities depend implic- 
itly on time via R(t) and 0(f). 

2.2 Torques and Inflow: Linear Regime 

When the perturbation is weak (\a\ <C 1), it is straightforward to 
expand the equations of motion for a test particle in the perturbed 
potential using the epicyclic approximation. With the substitutions 



R^Ro+Ri (R, t) =R +R a cos [m (ft - ft,,) t] 

-¥ 0o (0 + 01 {R, t) = o (f ) + (j> a cos [m (ft - Q p ) t] 

the linearized equations of motion become: 



(7) 



R { (R,t)=(^\ cos [m (ft - ft p ) t] 



fa(R,t) 

where 



-2ft 



i>n 



*i(0 

Ro R 2 (ft - ftp) 



cos [m (ft -ftp) t] (8) 



fd$ n 



+ 



2ft$„ 



dR fl(ft-ftp)J« 



A = K — in (ft — ftp) 



(9) 
(10) 



When the orbits of particles do not cross, i.e., there are no 
shocks, the net torque integrated over an orbit is identically zero to 
first order in $i . To see this, note that the instantaneous torque per 
unit mass r is 

9$ 



T = -(RX V$); = - 



d<p 



(11) 



Integrated over an orbit, the net angular momentum loss is therefore 



AJ- 



T&t 



9$ fl\-lAA 



(12) 



$->$o(fl) + $i(/f, 



(5) 



To first order, we must use = 0o = ft — ftp, so that AJ ~ 
-$ a (R) (ft - ftp) -1 f^{df{(t>)/d(f))d<j>, which trivially evaluates 
to zero because /(0) is periodic. 

Expanding the orbits to second order does not change this con- 
clusion - the net torque is again zero. Goldreich & Tremaine (1978 
1979 ) and Athanassoula ( 1992a , 2003 ) derive this result rigorously. 
They show that there can only be a non-zero second-order torque at 
the resonance points - where the epicyclic expansion breaks down 
- or if the perturbation magnitude is growing in time at a rate that 
is large compared to ft — ftp. 

The above derivation applies for collisionless particles. How- 
ever, gas is collisional, and thus obeys a slightly different disper- 
sion relation. For systems consisting of both gas and stars - which 
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are the focus of this paper - the gas mode is offset from the stel- 
lar mode by a small angle 8 due to dissipation. This gives rise to 
a torque on the gas by the stars (typically 8 ~ 10 — 15 deg in sim- 
ulations; |No guchi 1988] |Barnes~& Hernquist 1 996)|Barnes|1998| 
|Berentzen et al.|2007[ l~Following |Wada| ( |1994} , the phase shift 8 
and the resulting orbit-averaged torque on the gas can be derived in 
the limit of weak modes that do not induce orbit crossing or shocks 
in the gas: 



dt 



m 
~2** 



MR) 



v 2 



sin (8) B(R) 



where /„ is the stellar mass fraction in the disk, 



tan 



m k (fl p — fi) A 



B = 



(2[1 -n,,/nr' + din$ a /dinR) 2 

n- 2 ^/a 2 + m 2 A 2 K 2 (n-n. p y 



(13) 

(14) 
(15) 



and A ~ 0. 1 is a dimensionless constant that parameterizes the dis- 
sipation in the gas (formally, A is the dimensionless constant in a 
phenomenological 'friction' term in the momentum equation that 
is proportional to the radial velocity; the value ~ 0. 1 is motivated 
by comparison with hydrodynamic simulations). Note that the net 
torque is second-order in the perturbation amplitude \a\ ~ & a /V 2 . 
In addition, for a bar (m — 2) mode inside co-rotation, 8 > 0, and 
the gas mode leads the stellar mode (implying net angular momen- 
tum loss by the gas); this reverses outside co-rotation. 

Since the gas follows nearly circular orbits (7 ~ QR ), the 
orbit-averaged gas inflow rate follows trivially from the torque r = 
dJ/dt: 



M = 2nT, e 



,R 2 fl 



V 2 (l+dinV c /dinR) 



(16) 



Note the sign convention: positive M is outflow, negative is inflow. 
For the simple case of a constant- V c disk with &„ ~ \a\V 2 and typ- 
ical values for 8 ~ 10° calibrated from simulations, equation [T3] 
implies 



Mi 



amEgas-ft fi 



(l-Q p /n) 2 \2-m 2 {l-n p /Q) 2 



(17) 



where a ~ 2 — 3. Far from resonance, this reduces to M ~ 

M 2 M gas n. 

In addition to the torque between the gas and stars, non- 
axisymmetric waves can also directly transport angular momentum 
(with the coupling at resonances); this yields a qualitatively simi- 



lar torque that is also second-order in \a\ iLynden-Bell & Kalnajs 



|1972[ >. Kalnajs ( 1971) derives an exact solution for the logarithmic 
spiral, for which the inflow rate can be written 



M K 71 = - 



K 



d\K\ : 
d\kR 



n V3m , 
2^20 ' 



(18) 



where K = T[(m + ikR + l/2)/2]/T[(m + ikR + 3/2) /2] and k is 
the radial wavenumber of the mode. The resulting torque scales as 
\a\ 2 but is smaller than the torque between gas and stars (eg. |l3p 
by two factors: (1) a term ~ (c s /V c )Q~ ] ~ £ gas /£tot, and (2) a 
pre-factor d\K\~ / d\kR\ that depends on the structure of the mode, 
which scales as ^ \kR\/(i +m 3 ) for \kR\ « 1 or ~ \kR\~ 2 for 
\kR\ 3> 1. As a result, the predicted torques become small in both 
the global and local mode limits. 



2.3 Torques & Inflow: Non-Linear Regime (Orbit Crossing) 

2.3.1 General Criteria 

In the epicyclic approximation of the previous subsection, the per- 
turbation to the orbit of a particle is periodic in cj>. For a gaseous 
disk, however, this assumption breaks down if the perturbed orbits 
of particles cross each other, because orbit crossing in gas leads to 
dissipation (generally shocks) and inflow. Such dissipation breaks 
the periodicity in <f> that leads to the cancellation of the first order 
(in $ fl oc a) contribution to the torque in equation \\-\ - As a result, 
in the limit of strong orbit crossing an order of magnitude estimate 
of the inflow rate is 

Equation[T9]will derived more rigorously in what follows. 

Papaloizou & Pringle (1977) derived a useful criterion for 
when orbits cross. Consider two quasi-circular orbits at the same 
4>, with infinitesimal initial radial separation Sri. Using the lin- 
ear solution for each particle (eq. [8), the separation between 
the particles as a function of time is given by Sr(t) — Sri [1 + 
(dR a /dRo) cos(m^)]. Thus if 



\t\-\m inn i -I l ( dfo fo dA \\>-\ 



(20) 



then the orbits cross, i.e., 8r < 0, at some point (phase) in the orbit. 
We caution that this is a sufficient, but not necessary criterion for or- 
bit crossing. One can readily derive a similar criterion for when the 
particles cross in azimuth: j<90„/<90j > 1; and, clearly, if \R a \ > R 
there will be orbit crossing. Near-resonances, one can have orbit- 
trapping (where the gas librates about a specific radius), which will 
lead to shocks, for \dR„/dR\ as small as ~ 0.05 — 0.1 (see Binney 



& Tremain e|1987[ >. In general, even for linear modes, there is no 
purely local criterion for orbit crossing - one of the above condi- 
tions must be satisfied somewhere for orbit crossings to occur, but 
if so, it can in principle occur anywhere. A full census requires cal- 
culating the mode structure everywhere and the resulting particle 
orbits. Nevertheless we can make considerable analytic progress in 
several limits. 

2.3.2 Strong Orbit-Crossing Limit 

We now derive the net torque and inflow rate in the limit \C,\ 3> 1, 
in which the gas experiences strong orbit crossings and shocks. 
Prior to equation [20] we derived the pericentric separation of two 
perturbed particles at a given 4>. Generalizing this result, the sep- 
aration after moving from some arbitrary <pt to a final cf>f will be 
given by 8r(t) = <5r,[l + (dR a /dR ) (f(m<f> f ) - f(m<f>i))l where 
R[ — R a f(m<j)) represents the epicyclic expansion solution, i.e., 
f(m<f>) — cos(m<f)). Because £ 3> 1, there is orbit crossing at a 
range of <f>. We are thus not in a regime in which the epicyclic ap- 
proximation is valid in a global sense; however, we can consider it 
locally valid for a small amount of time ("between shocks"). Glob- 
ally, the gas orbits may be quite non-circular but we can approxi- 
mate this motion with a series of small epicyclic 'steps' separated 
by shocks that dissipate angular momentum and energy, allowing 
the gas to flow in. 

Given an initial separation of Sri at (/>,, the orbits will cross at 
a final <pf such that f(mcj>f) — f(m,(j)i) = —1/0 We assume that 
dissipation is efficient when the gas streams cross and shock be- 
cause the cooling time of the gas is always short under the condi- 
tions of interest. As a result, the colliding streams will dissipate the 
kinetic energy of their relative bulk motion. This effectively "re- 
sets" the orbit to a circular orbit with the appropriate total angular 
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momentum per unit mass of the streams at that time. We use this 
approximation to estimate the angular momentum loss due to the 
non-axisymmetric torques; the same results could be derived using 
energy considerations. 

The angular momentum per unit mass in the rest-frame is j = 
R 2 (4> + Q p ) = (Ro + R[ ) 2 (4>a + 4>i + Q p ) (recall that <j> is defined 
in the frame rotating at Q p ). Expanding to first order, the change in 
angular momentum of the streams from <f>j to <f>f is Ay = jj — y, ~ 
Ro (2R a Q +Ro<j}") (f[ m <t>f] — /[ m 0/]) where we have used R{ = 
R a f(ni(j)) and <f)\ — cj> a f(mcf>). Using equation Isl and the condition 
for when the streams shock (f[mcj>f] — f[m<f)j] = — 1/Q, we find 



A.r 



R[2R a ((f>0 + Qp)+Ro<pa} 



c 



*., 



1 



n-Q. p c 



The time required to move in azimuth by Acf> = 



Af = 



A(j> 
\Q-Q„ 



f-\f{m^)-C V )-m 



m\Q — Q D 



(21) 

(22) 
ij is given by 
^ (23) 



Note that because the f(mcj>) term ranges from — 1 to +1, the sign 
of £ = dRa/dRo in the above simply amounts to a phase offset in 
the shock location, so the above can be generalized by the replace- 
ment £ — > |C|. 

The average rate of angular momentum loss is simply given 
by Aj/At. Using the above solutions for At(<pt) and Ay, this can 
be written as a sum over the "number of shocks" encountered in a 
circuit of A(j> = 2-n. In the \C\ ^> 1 limit, we write this as continuous 
average over cf>, with Aj/At — > (dj/dt) — r and 

Ay 



<r) = 



-d4>i 



(24) 



At[<f>i] 

where (r) is the azimuthally averaged torque. The general expres- 
sion for t cannot be written in closed form, but the leading-order 
term in £~ can be used to derive the following expressions for the 
torque r, inflow rate M, and radial velocity Vr\ 



(r) 



M - 



2tt 
(E, 



$„F(C)sign(fi-n„ 



S .R 2 S1)[ 



msign(S7 — Q p 



V R 



H0> 



*[§ 



(l+d\nV c /dlnR' 
msiga(Q-Q p ) ^ 



no 



(l+dlnV c /d\nR) 



2- 



2|C 



+o(\cr 2 ) 



sign(C) (|Cl»l) 



(25) 
(26) 
(27) 
(28) 



Note that the net torque and inflow rate are independent of £ 
for ( 3> 1. They depend primarily on the magnitude of the non- 
axisymmetric potential <!?„. This is because although the change in 
angular momentum per shock scales as £ -1 (eq.|22[>, the time be- 
tween shocks itself also scales as £ _1 (this follows from Taylor 
expanding eq.[23]and averaging over </>,)■ Another way to say this 
is that there are N ~ ( shocks each of which resets the gas to a cir- 
cular orbit with Ay oc 1/f, so that the net torque is independent of 

c 

Dimensionally, given <J>„ < 



26 



reduces to M > 



\a\ V c , equation^ 
— \a\ E gas ^ 2 Q, in agreement with our order of magnitude estimate 
in equation [19] Thus in the strong-torque regime the inflow rates 
are linear in \a\: when orbit crossings are important, gas inflow can 
easily be enhanced by an order of magnitude or more relative to the 
expectations from the weak-torque theory of §2.2| 

The sign() terms in equations |25|28| are important and 
arise from terms that in the full solution are continuous; i.e., 



G/\/l + C 2 ~> sl g n (C) f° r ICI 3 s 1- Recall that with our sign con- 
vention, $„ is negative where E„ is positive. For typical parameter 
values inside of co-rotation (SI > Q, p ), £ is negative and therefore M 
is also negative (inflow), whereas outside of corotation it is positive 
(outflow). This dependence on corotation is essentially a question 
of whether the shocks are in prograde or retrograde motion in the 
pattern frame, which changes at the corotation resonance. The de- 
pendence on sign(£) is due to the phase between the shocks and 
the surface density perturbation: when £ < 0, shocks occur where 
f(m(j)) — f(m<f>i) is maximized - in phase with the positive max- 
imum in Ei - and hence where (for Q > Cl p ) Ay is positive. This 
generically leads to outflow. But when £ < 0, the shocks occur lead- 
ing the mode, where Ay is negative and minimized with respect to 
the unperturbed orbit; this leads to inflow. 

If instead of using real variables, we generalize to complex <!>„ 
and (, a nearly identical derivation to that presented here leads to 
<I> ( ,sign(^) — > |$ | cos(arg[$ C*]), i.e., |$ a | cos0, where is the 
phase angle between £ and $ a ; this can have an in-phase inflow 
and/or out-of-phase outflow term. Likewise, the sign(f2 — fl p ) term 
can be generalized for complex u> = mil p + iy, by taking sign(fi — 
Qp) — > Re[m\Q — Q p \ / (mQ — to)]. 

2.3.3 Marginal Orbit-Crossing Limit 



The derivation in i 2.3.2 is only valid when \Q 3> 1. Indeed, when 
ICI < 1, the quantities used in these derivation are undefined, even 
though it is still very much possible that there are shocks and even 
orbit crossings (as discussed in § |2.3. 1 1 >. Moreover, the requirement 
of formal orbit crossings for strong torques and shocks is a conse- 
quence of the fact that we have thus far (implicitly) considered only 
infinitely cold disks. There can, in fact, be shocks when \(\ < 1, 
even without orbit crossings, if the non-axisymmetric potential is 
dominated by the stars rather than the gas. We now estimate the 
resulting torques and inflow rate in this limit. 

If the disk is not gas-dominated, then the gas orbits are given 
roughly by the standard equations of motion in response to the stel- 
lar potential. Given radial variations in some property of the gas, 
the length-scale /z e ff at which p ressure forces become important is 



given by A c ff = 2c cff /GE (e.g., Toomre 19641 where c c ff is the ef- 



fective dispersion of the gas, including both the true thermal sound 
speed and any non-thermal turbulent motions: c 2 ff — c 2 + v 2 . Note 
that in a spherical potential, h e ff is different from the vertical scale 
height h by ~ 2c e ff Q/GT, ~ Q, the Toomre Q parameter (defined 
using c c ff not c s ). The stellar motions do not explicitly appear in the 
effective sound speed here, because the key criteria is whether gas 
shocks; E does, however, represent the mass density of gas plus 
disk stars since both participate in the mode self-gravity. 

If stellar streamlines approach to within a distance < hes, they 
remain collisionless and conserve energy. But if gas that is ini- 
tially spatially separated is forced into a supersonic 'near-collision' 
within a distance < /i e ff by the potential of the stars, it will shock 
and radiate the difference in energy between the streamlines (con- 
serving only their total momentum). The key point is that when the 
non-axisymmetric potential is stellar dominated, gas "collisions" 
can happen not just at strict orbit crossings, but at a finite separa- 
tion ~ /j e ff between streamlines. 

If two gas streamlines approach one another at a relative ve- 
locity that is much smaller than the sound speed c s , pressure effects 
can readjust the orbits without much dissipation. Thus significant 
dissipation probably requires relative stream velocities > c s . What 
are the relative velocities? For radial motions, the maximum mag- 
nitude of the relative velocity over an orbit, from the epicyclic so- 
lution, is Av = \dvijtfdR\8n = [h^/{l - \Q)]\dR\/dR\, where 
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in the second expression we have used the maximum initial sepa- 
ration streamlines can have such that they are separated by /i e tf at 
some point in the orbit, Sri = /j c ff/(l — |CD- The resulting criterion 
for gas streams to approach within a distance ~ h s a with supersonic 
relative velocities Av > c, is given by 



1 dvi» 



OR 



Ceff 
Cs 



1 - ICI 

2-kQ ' 



(29) 



Since vi,» is linear in the mode amplitude \a\, we can define 
\dvi.R/dR\ a =i as the magnitude of dvi,ji/dR for \a\ = 1, and write 
equation[29]as a condition on the mode amplitude 



H> 



(i-lCI ) (c 

2nQ 



Cell 



r ■>( 


dvi.R 


2 


dv\.<i> 


2 w 


n 


dR 


+ 


cm 




L \ 





-1/2 



(30) 

where we have also included the relative azimuthal velocity vi,^ 
for completeness. 

In equations [29] and [30] we have allowed for the possibility 
that the disk has non-thermal random velocities via the effective 
sound speed Ceff, which is also in our definition of Q. One sub- 
tlety is whether the condition for effectively dissipating the non- 
axisymmetric motions is a relative velocity Av > c e s or Av > c s . 
We suspect that the latter is appropriate. When CeS > Av > c s , the 
energy per unit mass dissipated via non-axisymmetric motions is 
small compared to the overall rate at which supersonic turbulence 
is dissipating in the ISM (the latter being due to, e.g., star for- 
mation). Nonetheless, the ordered velocities induced by the non- 
axisymmetric potential can be efficiently dissipated if Av > c s \ 
we confirm this below in simulations that have resolved turbulent 
Ceff > C s . 

Equations |29| and |30| cannot be satisfied for \a\ < 1 in a suffi- 
ciently cold disk with 2< 1 and |C| < 1. In this limit, the condi- 
tion for shocks and orbit crossing reduces to our previous criterion, 
C > 1 . However, for larger Q, the effective "thickness" of a given 
gas stream increases and so for the same non-axisymmetric pertur- 
bation, it is easier to generate collisions of streamlines. Of course, 
for sufficiently large Q our assumption that the gas orbits are well 
described by test particles responding to the stellar potential will 
break down. But for the particularly relevant case of Q > 1 and 
t'eff > c s , equations [29] and [30] show that shocks induced in the gas 
by the stellar potential are possible even though |C| < 1. For typ- 
ical modes (say, m = 1 modes in Keplerian potentials or bars in 
constant- V c disks), at the order of magnitude level, this condition 
on the mode amplitude becomes \a\ > 0.2(2 ' (c s /Cea). We reit- 
erate that this is only true when the gravity of the perturbation is 
dominated by stars. For gas dominated systems, the mode structure 
would already include pressure effects, invalidating our analysis. 
The exact gas fraction at which this transition occurs is sensitive 
to the gas equation of state and phase structure; it is thus diffi- 
cult to estimate precisely, beyond the order of magnitude estimate 
M g - ds >M t . 

We now generalize the derivation in j p.3.2| to determine the 
inflow rate and torques on the gas in the |C| < 1 limit. We as- 
sume that both C and Jua/R are small, so that we can Taylor ex- 
pand in both quantities. The separation between two streamlines 
is given by 8r((f>) — Sri (1 + C [f(m<f>) — /(m0,)]) Gas that is ini- 
tially within a maximum separation <5r, = SR mdx = hes/(l — |C|) ~ 
Aetf(l + ICI) wl ll collide at some point in the orbit. However, stream- 
lines with initial separations less than a minimum value SR mm 
are never differentially accelerated by the perturbed potential to 
Av ~ c s and thus never experience shocks. The exact value of 
8R m [ n depends on |a|, (, |ch>i,*/9/?| a =i, etc. As a rough estimate, 



we take SR mm ~ /? e ff because gas initially separated by SR m [ n <■ /feff 
is strongly affected by internal pressure forces and is unlikely to 
be accelerated to relative velocities > c s . For a given <5r, satisfying 
SRmin < Sri < 5R miix , the streamlines reach Sr < h^g (i.e., 'collide') 
when f(m<f>) — f{m<f>i)-\-C > ~ 1 {h e s/5ri — 1). For<5r, =^ m in, this cor- 
responds to f(m(f>) = f(m<j)i) and A<f> = 0, while for <5r, = R mdx , 
the shocks occur at f(m<f>) = f{m(f>i) ± 1 and A(f> — ±(ji/2m). 
When the gas reaches the <j> such that Sr(4>) < fhs, it dissipates 
its relative orbital energy but conserves momentum. From j |2.3.2[ 
the change in specific angular momentum Ay for streamlines satis- 
fying SRmin < Sri < SRmta is given by 



Aj- 



$ fl 1 



n - fi p c 



fleff 

~SJi 



(31) 



The average change in angular momentum over all initial stream- 
lines separations is given by 



<Ay) : 



1 



frefflCI 



Wl + ICI) 



Ajd(Sn) 



1 $ a 



2 n-n, 



sign(C) (32) 



where we have used R m [ n = h e ff and R mdK = fr e ff(l + |CI) an d in 
the second equality have Taylor expanded for C <C 1 . Note that /i c ff 
has cancelled out in our result for (Ay) - the value of the stream 
thickness does not influence how much angular momentum is lost. 
The integral over Sri in equation[32]is equivalent to an integral 
over the cf> of the shocks. Thus the average over Sn is equivalent to 
an average over cj> and/or time. In particular, the range Sri = R m m 
to Ram. corresponds to f(mcj>) — to ±1. This is equivalent to an 
integral over time for a stream circuit of (j> — 2iv/m; which takes a 
time At — 2n jm \fl — fl p \ . To leading order, we therefore obtain 



<r) = — <f- fl F(C)sign(r2-f2 p ) 

Z7T 



M=C£ g ^R 2 n)[^ 



$ a ] msign(fi — Q p ) 



{l + dlnVc/d\nR) 



no 



/■-iu^Mgnio^+oacr 1 : 



(33) 
(34) 

(35) 



Note that equations |33|35| in the low C limit are comparable to equa- 
tions |25|28] derived in the C 3> 1 limit (to a factor of ~ 2), subject to 
the additional restriction that equation [30] for the mode amplitude 
\a\ (and thus C) must be satisfied in order for there to be shocks in 
the low C limit. 

2.3.4 Combined Scaling 

For |C| ~ 1, we can reasonably interpolate between the two regimes 
considered in the previous two subsections (|C| «C 1 and |C| 3> 1) 
to find: 



M=(E gas ^ 2 r2) 

S(w,$i)=Re 

1 



Ijal m5(cj,$i)F(Q 
IK 2 ! l+dhiVc/dinR 



m | Q. — Q, p 
mil — uj 



ICI 



F(C) ^2 + 6 + 2|C|/3 



+ - 



cos(arg[$„C*]) 

ICI 2 



(36) 



(37) 



(38) 



1 + (|C|+9|C| 2 + |CI 3 )/13 
Recall again that in the low C limit, these expressions require that 



equation 30 for the mode amplitude \a\ be satisfied. 



Note some critical features of these results. The inflow rate 
M is linear in |a| and thus, as advertised, shocks induced in the 
gas by the non-axisymmetric stellar potential (when present) will 
dominate the total angular momentum loss relative to the direct 
transport by the gaseous mode itself. We collect the phase infor- 
mation in the term 5 in equation [37] (5 < giving inflow, 5 > 
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outflow). If the modes are purely real and global S — > sign(fi 



Clp) sign($ a ) sign(£). The generalization of 5 in equation 37 allows 



for modes with radially varying phases, i.e. <!>„ oc exp (i J kdr), as 
in j ]2.3.2| It is easy to show that for a mode in a galaxy with typical 
structural properties, there is generally inflow induced inside of co- 
rotation, and outflow outside (for sign($ a ) < 0, then sign(f ) > 
and sign(fi — Q p ) changes from positive to negative from inside 
to outside corotation); this is consistent with the shock-free, weak- 
torque results in |2.2| Also for typical structural parameters, the 
phase changes at the Lindblad resonances, producing outflow in- 
side of the ILR and inflow outside of the OLR. 



3 EQUILIBRIUM GAS DENSITIES 

In 5J2] we described the physics that sets the inflow rate of gas due 
to gravitational torques. Statistically, the net inflow rate through 
a given radius in a galaxy is set by a competition between these 
torques and star formation. This competition determines the net in- 
flow rate at very small radii (e.g., onto a central BH) and the result- 
ing gas surface density profile. 

According to the observed Schmidt-Kennicutt law, the average 
star formation rate surface density is given by 



E* 

tK 



(B 



(39) 



where the constants t)k ~ 3/2 — 7/4, Ek and tK are discussed 
below. In the limit of a very small supply of gas, the non- 
axisymmetries in the potential will in general be modest and the 
inflow of gas will not be able to compete with star formation, which 
will consume most of the gas. In the opposite limit of a large supply 
of gas, however, the inflow rate approaches a linear function of E gas , 
while the star formation rate is a nonlinear function. Thus the gas 
will be consumed rapidly until the inflow and star formation reach 
an approximate, time-averaged equilibrium. We will consider this 
in more detail in j ]5.2| but for now note that inflow and star forma- 
tion balance each other at a quasi-steady gas density (E eas ) given 

by 



<E gas > 



\M\tK \Vi* 



nR 1 T. K 



(40) 



Moreover, because M itself depends on E gas , equation 1 40 1 implies 



an equilibrium gas surface density profile. For the simple dimen- 
sional scaling \M\ ~ \a\ T} gim R 2 Q., we find 



E gas (i?)~E*(^fi(.R) 

TV 



i/Ow-i) 



(41) 



Since fi ~ \jGM mc (< R)/R 3 , this implies that the scaling of the 
total surface density sets the equilibrium gas density profile, if we 
are in the strong torque regime. 

We stress that equation [41] and the associated balance be- 
tween star formation and inflow are only approximately valid, and 
only then in a time-averaged, order-of-magnitude sense. There are 
regimes (e.g. the early stages of coalescence in a major merger) 
where gas can flow in on of order a single dynamical time (Barnes 
& Hernquist 1991 , Mayer e t al.|2010) , much faster than star forma- 
tion (observationally supported by observations of e.g. ULIRGs 

[Soioirion et al.|[T997[ 



with very high nuclear gas fractions; see 



Brya nt&Scoville|1999||Tacconi et al.|1999| 



Hibbard&Yunl999). 



In the opposite extreme, if gas densities are sufficiently low, the 
star formation and/or inflow timescales can be much larger than the 
dynamical time. We show below that simulation gas densities typ- 
ically scatter by an order of magnitude about the result implied by 



equation[4T] Nonetheless, this is a useful estimate of the gas densi- 
ties achieved, particularly for analytic purposes. It is also important 
to note that the equilibrium between star formation and inflow still 
allows for an arbitrary range in observed gas fractions, as new stars 
are continuously produced. 

To use equation [41] for specific predictions, we require val- 
ues for E K , t K , and r\ K . From the fit in Kennicutt (1998), t K ~ 
0.63 x 10" yr and r) K « 1.4-1.5, for E* = lO 8 M kpc -2 . How- 
ever, th ere is considerable uncertainty. For example, Bouche et al. 



( 2007 i suggest t K m 0.41 x 10*yr and r/ K « 1.71 ±0.05, for E* 



10 MQkpc~", from observations of high-redshift systems and 
more extreme starbursting galaxies. Davies et al. (2007 ) and Hicks 
|et al.| p009 ) measure star formation in nuclear regions around AGN 
(~ 1 — lOOpc), perhaps closest to the regions of interest here, and 
find results consistent with the scalings in Bouche et al. (2007). 
When necessary, we use the latter throughout this paper. 

4 COMPARISON TO HYDRODYNAMIC SIMULATIONS 
4.1 The Simulations 

Thus far, we have made a number of simplifying assumptions in 
order to derive tractable analytic results. In this section, we com- 
pare these analytic models to numerical simulations of gas inflow 
in galaxies. We do not present any new simulations, but instead 
draw on results that have been described previously in a number 
of papers (see e.g. Di Matteo et al.|2003 Robertson et al.|2006b; 



|Cox et al.|20"06]|Younger et al.|2 008. Hopk ins & Quataert|2 010a). 
The simulations we compare to include galaxy-galaxy mergers, iso- 
lated galaxy disks, and multi-scale "re-simulations" that follow gas 
from large scales in galaxies to sub-pc scales around a BH. TablefT] 
summarizes the broad properties of these simulations and provides 
references to the original papers. 

The simulations were performed with the parallel TreeSPH 
code GADGET-3 ([Springel 2005}, based on a formulation of 
smoothed particle hydrodynamics (SPH) that conserves energy 
and entropy simultaneously even when smoothing lengths evolve 
adaptively (see e.g., Springel & Hernquist 2002; Hernquis t|1993| 
O'Shea et al. 20031. The detailed numerical methodology is de- 
scribed in iSpringel (2005), Springel & Hernquist (20031, and 
Springel et al. ( 2005 . The model galaxies consist of gas disks, stel- 
lar disks, stellar bulges, BH particles, and dark matter halos, with 
Hernquist (1990) halos, exponential disks, and (optional) bulges 
with mass fractions B/T ~ 0.1 — 1. Some of the simulations in- 
clude sub-resolution prescriptions for BH accretion and feedback, 
but they are not important for the physics described here - for our 
results, only the BH mass is critical, and then only on the smallest 
scales < 10 pc where it dominates the gravitational potential. On 
all scales, extensive resolution tests have been performed (see ref- 
erences in TablefT), with high-resolution case studies using > 100 
times as many particles. 

We consider a large number (~ 50) of simulations of isolated 
galactic disks, marginally to strongly self-gravitating (described in 
detail in Younger et al. 2008). We also consider a large number of 
mergers of these idealized disks, in which we vary the merger mass 
ratio (from 1:20 to 1:1), orbital parameters, and structural prop- 
erties of the disks (for details, see Hopkins et al. 2009b). These 
structural variations include changing the disks in mass, bulge-to- 
disk ratios, disk scale height, gas and stellar disk scale lengths and 
mass profiles, and scaling the disks following Mo et al.] p998) to 
correspond to expected disk properties at various redshifts from 
z = — 6. Finally, we also include the "re-simulations" of galactic 
nuclei described in detail in Hopkins & Quataert (2010a). In these 
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Table 1. Overview of the Simulations 



Simulations 


Spatial Scale 


Mass Scale [Mq ] 


e 


[pc] 


Major Mergers 


0.1- 


- lOOkpc 


10 lJ - 


-10 12 


10 


-100 


Minor Mergers 
Isolated Disks 

Intermediate-scale Re-simulations 
Nuclear-scale Re-simulations 


0.1- 

0.1 

10- 

0.1 


- lOOkpc 

- lOkpc 
lOOOpc 

- lOOpc 


10 9 - 

10 9 - 
10 7 - 
10 6 - 


-10 12 
-10 12 
-10 10 
-10 8 


10 


-100 
-30 
-1 
0.1 
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Younger et al. 1 2008 1, Tables 1-2 
Younger et al. 120081 Tables 3-4 
Hopkins & Quataert (2010a , Tab 



(2006a 



, Table 1 



Hopkins & Quataert (2010a , Table 



Cox et al. 



es 1-2 

3 



(2008 



List of SPH simulations used to test our analytic predictions. From left to right, columns summarize: (1) Simulation "class," each of which includes 
~ 100 individual simulations, with different structural properties (mass profiles, kinematics) of the initial disks, bulge-to-disk ratios, BH mass, gas 
fraction, baryonic mass, orbital parameters (in mergers), and efficiency of stellar feedback and star formation. Particle numbers range from ~ 10 5 — 10 7 
(2) Approximate range of spatial scales spanned by each simulation class. (3) Corresponding total range in baryonic mass. The particle mass is ~ 10 — 6 
times smaller. (4) Typical spatial resolution/softening length (often varied with galaxy scale length). (5) Reference(s) for the details of each simulation 
class. Particular tables and figures in the original references that summarize the simulation properties are noted. 



calculations, we re-simulate the central ~kpc of a representative 
sub-set of the mergers and isolated disks at several times in their 
evolution, at much higher resolution, for a total of ~ 100 simula- 
tions. This technique allows us to study, using realistic initial condi- 
tions, the formation of secondary instabilities when the gas/stars are 
self-gravitating ("bars within bars") and the resulting gas inflow to 
small radii. These simulations are accompanied by a similar num- 
ber of yet smaller-scale re-simulations that follow the gas inflow 
deep into the BH potential, with ~ 0.1 pc resolution (the resolu- 
tion studies mentioned above include both varying this spatial res- 
olution and the scales at which the re-simulations are performed; 
neither of these changes our conclusions). 

In both sets of "re-simulations," we consider simulations with 
zero initial perturbations, with initial perturbations seeded ran- 
domly, and with initial perturbations inherited from the "parent" 
simulation. In Appendix A3 of Hopkins & Quataert (2010al we 
show that this does not qualitatively change our results, because 
the systems of interest are self-gravitating on small scales and so 
the local structure determines the equilibrium modes, not the ini- 
tial conditions. Moreover our general analytic model appears to de- 
scribe shocks due to modes from either initial condition equally 
well. 

The simulations utilized here all include gas cooling and star 
formation, with gas forming stars at a rate motivated by the ob- 
served Kennicutt ( 19981 relation. Specifically, we use a star forma- 
tion rate per unit volume p\ oc p 3 ' 2 with the normalization cho- 
sen so that an isolated Milky-way like galaxy has a total star for- 
mation rate of about lMQyr -1 . This formulation produces a pro- 
jected Schmidt-Kenicutt relation (i.e., E*) similar to that measured 
by Bouche et al. (2007). We crudely account for feedback from 
stars with an effective equation of state that mimics turbulence in 
a multi-phase interstellar medium by adopting a sound speed that 
is much larger than the true thermal sound speed (e.g., Springel & 
Hernquist 2003 1. Hopkins & Quataert ( 20 1 0a ) describe in detail the 
effects of different subgrid ISM sound speeds on angular momen- 
tum transport and inflow rates. They also argue for effective turbu- 
lent speeds of ~ 10 — 50 kms -1 for densities ~ 1 — 10 5 cm -3 , re- 
spectively. As we have shown in €2] the torques and inflow rate in- 
duced in the gas by an asymmetric stellar potential are independent 
of the ISM sound speed in the orbit-crossing limit (e.g., eq. |36[ l. 
This suggests that the use of a subgrid sound speed in the simula- 
tions does not have a significant effect on the physics of gas inflow. 
However, the critical assumption in both the numerical simulations 
and the analytic calculations presented here is that stellar feedback 



maintains a reasonable fraction of the ISM in a diffuse phase, rather 
than being bound in star clusters. In the simulations, gravitational 
instability can form some compact, dense clumps (analogous to 
molecular clouds or star clusters), but the enhanced subgrid sound 
speed ensures that most of the ISM does not suffer this fate. This 
is intended to mimic physics not included in our calculations, such 
as the observed low efficiency of star formation in dense molecular 
gas (e.g., Krumholz & Tan 2007) and the disruption of molecular 
clouds by radiation pressure, HII regions, etc. (e.g., Murray et al. 
|2010[ > - the latter returns most of the gas from molecular clouds 
back into the diffuse ISM. 

There are of course considerable uncertainties in any model of 
the ISM and star formation in galaxies, especially on small scales 
near a BH; the adoption of a simple Kennicutt (1998) relation is 
merely a convenient simplification. For this reason, |Hopkins &| 
Quataert (2010al consider a broad range of such prescriptions (see 
their § 2 and Appendices). We briefly summarize some of these 
variations, but interested readers should see that paper for details. 
The slope of the star formation law was varied from p» oc p'~ 2 , 
amounting to changes in the star formation efficiency on small 
scales by factors of > 100. This spans the entire observed parame- 
ter space of star formation efficiencies on large galactic scales (see 
e.g. Bigiel et al. 2008 ), in dense molecular clouds ( Krumholz & Tan 
|2007| >, and the central regions of AGN l |Hicks et al.|2009| >; our "typ- 
ical" cases are consistent with the latter observations. The primary 
difference the star formation efficiency makes is that it changes the 
gas-to-stars ratios at a given time; but since this quantity explicitly 
appears in our model, the effect is to move systems along our pre- 
dicted scalings, not with respect to them. In our analytic estimates 
that follow, when the absolute value of star formation efficiency is 
important, we parameterize it as a free function, allowing for dif- 
ferent detailed models of star formation. 

The simulations in TablefTlalso include large systematic stud- 
ies of the prescription for the sub-grid gas physics in the ISM. 
Specifically, the effective equation of state of the gas was widely 
varied. The simulations include cases where it is quite stiff on small 
scales (near-adiabatic, similar to what is adopted in the studies of 
|Mayer et al.|2007||Dotti et al.|2009) , through to cases where the gas 
is allowed to cool to a cold isothermal floor (with c s <C V c ), and thus 
forms a clumpy, inhomogeneous medium (on galactic scales, simi- 
lar to what is assumed in Bournaud et al. 2007 ; Teyssier et al. 2010). 
Intermediate cases are also sampled, allowing for e.g. an equation 
of state similar to that presented in Klessen et al. ( 2007 1, motivated 
by the models of Spaans & Silk (2005 ) for dense star-forming gas. 
In a large number of our nuclear-scale simulations, the gas disper- 
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Figure 1. Inflow properties in a typical gas-rich (/ gas « 0.4) merger of two equal-mass ~ Z,* (M 4 = 10"Mq) galaxies (run b3ex in Hopkins & Quataert 
2010a), at different times: before any interaction, when the system is unperturbed (left); just after an early first passage (center left), which induces a bar 
response in the primary; somewhat later, just before second passage, as the induced bar feature has relaxed (center right); just after final coalescence (right). 
Top: Images of the gas density at each time. Intensity reflects the gas surface density. Color reflects the specific star formation rate (from blue being low, with 
gas consumption timescale 2>Gyr to yellow being high, with consumption timescale ~ 10 7 — 10 s yr). Horizontal axis is cartesian (x, in kpc). Second from 
Top: Gas and stellar surface density profiles as a function of cylindrical radius R. We also show the analytic gas density profile from eq.|4T] which assumes that 
gas inflow approximately balances star formation at each radius (see §151. Middle: Amplitude of non-axisymmetric modes versus radius. Second from Bottom: 
Orbit-crossing measure |C(^)| = \3R\/dR\ (eq |20} ; \C\ > 1 implies strict orbit crossing, but shocks are still present at many radii with |£| < 1. Bottom: Gas 
inflow (negative) or outflow (positive) rate as a function of radius. Dotted red line shows M glis = 0. We compare the simulation results to our analytic prediction 
based on strong shocks induced in the gas by t he n on-axisymmetric stellar potential (§ |2.3.4) . The predictions from the no-shock, weak-torque limit (eq |17) or 
the direct transport by the mode in the gas (eq|l8i are smaller by a factor > 10, and are would be indistinguishable from the M = line. 



sion is dominated by resolved turbulent gas motions, which main- 
tain an effective c s even as the sub-grid c s becomes very small (the 
origin of this turbulence and its importance for e.g. AGN torii will 
be discussed in future work). As such, we can explicitly check that 
our formulation for gravitational torques still applies even in turbu- 
lent inhomogeneous gas distributions. 

Nevertheless, we consider our sub-grid treatment of the ISM 
to be the biggest uncertainty in the numerical models. Because we 
cannot resolve the extremely large densities of actual star-forming 
cores in these regions, star formation is almost certainly more dif- 
fuse/uniform throughout the disk, which may in turn have conse- 



quences for the ability of stars and gas to interact efficiently. And 
the ISM is, by definition, more smooth than the implicit sub-grid 
multiphase structure. An important goal for future work will be the 
investigation of more realistic ISM and star formation models on 
these scales. These are unlikely to make a large difference in the 
limit of stellar-dominated disks (our focus here), but may in the 
very gas-rich regime open up new channels for angular momentum 
exchange. 



In 34721 we provide an overview of the simulation results and 
a comparison to the analytic results derived in ip] & [3] for exam- 
ple simulations covering three different spatial scales: galactic (~ 
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Figure 2. Analogous results to FigurefTlfor re-simulations of the central ~kpc in galaxy mergers, at representative times in the large-scale inflow history (from 
left to right, these are runs If9b5res, Iflbllate, Ilowresq, Inf28b4 in Hopkins & Quataert 2010a). These simulations resolve gas inflows to within ~ lOpc. 



kpc), intermediate-scale (~ 0.01 — 1 kpc), and nuclear (~ 0.1 — 10 
pc). We refer the reader to Hopkins & Quataert (2010a) for more 
detailed information about the properties of the gas and stars in the 
simulations (e.g., Q[r], h/r, etc). Although the specific simulations 
discussed in j |4.2| are quite instructive, it can also be misleading 
to focus on the results of an individual simulation. The reason is 
that there is large variation in time and potentially large scatter in- 
troduced by modest differences in galaxy properties (due to, e.g., 
large-scale fragmentation of a spiral arm biasing the results of a 
given simulation). For this reason we believe that it is important to 
consider the statistical properties of a large number of simulations. 
In § |4.3| we thus present a more comprehensive statistical compari- 
son between the simulations and our analytic results. 



4.2 Overview of Simulation Results 

Figure [T] provides an overview of the dynamics during a major 
merger of equal-mass ~ L» galaxies with / gas « 0.4. The first row 
shows, from left to right, images of one of the two galaxies before 
first passage, after a near-final passage (when the impact parameter 
of the interaction is < the size of the galaxy), and just before the 
final coalescence of the two systems; the final column shows an im- 
age of the remnant system just after coalescence. The second row 
in Figurefllshows gaseous and stellar mass surface density profiles 
at each time, while the third shows the amplitude of the stellar sur- 
face density perturbations for different azimuthal mode number m 



as a function of radius (R — is defined as the location of the BH); 



t (R,t)\ = 



(42) 



|/ 27r E(i?,<£)exp(fm</>)d<A| 

The phase angle 4>q of the mode maximum is determined using 
a similar expression, as are the potential perturbations, and other 
perturbed quantities of interest. The second to last row in Figure [T] 
shows estimates of the orbit crossing parameter Q for the m = 1 
and m — 2 modes (the calculation of which is described below). 
Finally, the bottom row in Figure [T] shows the instantaneous gas 
inflow rate as a function of radius R, determined directly as M — 
AR~ J dMgas vr given the radial velocity vr of each gas particle of 
mass dMg as in a small radial annulus AR. As before, we define M < 
to be accretion/inflow, whereas M > corresponds to outflow. 

In addition to the numerical results, Figure [T] also shows sev- 
eral of the predictions of our analytic model. First, given the mea- 
sured inflow rate in the simulations M, we use the balance between 
inflow and star formation implied by equationHOlto predict the gas 
surface density E gas (2nd row). This prediction is within a factor of 
~ 3 of the measured gas density in the simulations, particularly at 
later times (right-most two columns) when our assumed balance 
between inflow and star formation is a better approximation. In 
addition, given the measured \a\ at each radius, and the local gas 
properties (E gas , fl, and their derivatives), Figure [T] also shows our 
prediction of the inflow rate M usin g eq uation |36| (bottom row). 
Evaluating the sign of M in equation |36J (i.e., inflow or outflow), 
is non- trivial given its dependence on Q p and f; one option would 
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Figure 3. Analogous results to FigurefTlfor re-simulations of the central ~ 10 — 50pc in the intermediate-scale simulations shown in Figurep] at times near 
the peak of accretion and gas inflow (from left to right, these are runs Nf8hlclqs, Nf8h2b4, Nf'8hlc0hol, Nf8hlcldens in Hopkins & Quataert 2010ai. These 
simulations have a resolution ~ 0. 1 pc and extend deep into the potential of the central BH. The (non-wound) m = 1 eccentric disk and/or single-armed spiral 
mode is ubiquitous in the quasi-Keplerian potential. 



be to fit a global model to the mode structure as a function of time 
and radius, but since £ depends on second derivatives this is very 
noisy. Instead, we found it easier and more numerically robust to 
use the epicyclic approximation to estimate £ and the sign of M. 
With this approximation, the mean local radial velocity (zero in the 
unperturbed state) obeys 

dvi.R 



dR 



■ im(n-n p )c, 



(43) 



i.e. it is phase-shifted by n/2 and multiplied by (fi — Q, p ) rel- 
ative to f; this implies that sign(£7 — Q,,) cos(arg[$i £*]) = 
— sin(arg[$i] — arg[dvi/dK\). We use the true radial velocity of 
the stars in the simulations vr as our estimate of the perturbed ve- 
locity vt ,r and determine its amplitude and phase as in equation [42] 
Together with the measured $i in the simulations we can thus de- 
termine the sign of M induced by a given m. For the results in Fig- 
ure [jj we sum the contributions to M from m — 1 to m = 5, but find 
that this is almost always dominated by the largest mode or two 
(this will be an assumption of ours later). Equation[43]also implies 
that \(\ ~ \dv[,R/dR\/[m\Cl— Q p \]; we estimate Q p from a global 
fit to the mode amplitude and axis cf>o as a function of time, assum- 
ing flp is constant in time and radius for a given simulation. The 
results for |£| are shown in the 2nd to last row in FigurefTlfor the 
major merger simulations. 

Figure [TJ shows that the predicted inflow rate as a function 
of radius is well reproduced by our analytic model at all times, 



providing support for the simplifying assumptions we made in the 
analytic derivation. By contrast, the predictions of the no-shock, 
weak-torque limit (eq|17[) or the direct transport by the mode in the 
gas (eq |18) are smaller by a factor > 10. Although the results in 
Figure [JJare for a major merger, we find similar consistency be- 
tween the simulation results and our analytic predictions for both 
minor mergers and isolated, bar-unstable galactic disks. Figure [TJ 
also shows that the orbit crossing parameter £ is ~ 1 at a range of 
radii, with |£| ~ 10— 100 3> 1 near the co-rotation resonance, im- 
plying strong orbit crossings and shocks. As described in the text 
above equation[43] our estimate of £ is itself based on the epicyclic 
approximation; it is thus very approximate in mergers given that 
the galaxies are strongly disturbed and quantities such as Q,, are 
not that well-defined. Nonetheless, the results in Figure[TJare con- 
sistent with the presence of strong shocks and dissipation in the 
simulations, which drive the large gas inflow rates. 

Figures [2] and [3] show analogous results to Figure [TJ but for 
simulations of the central ~kpc (Fig. [2| and ~ 10 pc (Fig. [3l of 
gas-rich mergers (from Hopkins & Quataert 2010a). The simula- 
tions span different morphologies and degrees of dumpiness in the 
ISM gas, and include both cases that are strongly stellar-dominated 
and that have more comparable gas and stellar densities. In each 
case, the agreement between the analytically predicted E gas and 
M and the numerical results is again quite good. As noted above, 
we do not necessarily expect our predictions to be applicable in 
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Figure 4. Effects of simulation resolution on the modes we compare to our 
analytic model. We plot the m = 2 mode amplitude (top), |C(^)| f° r the 
m = 2 mode (middle), and inflow/outflow rate (bottom), for intermediate- 
scale simulations as in Figure[2] For each, we show three otherwise identical 
simulations: one with low resolution (/Vpanicles = 5 X 10 5 ; comparable to the 
resolution in the central regions of our galaxy-scale runs before "zoom-in"), 
one with our standard resolution (N = 5 X 10 6 ; that of the other zoom-in 
runs in Figure |2|, and one with very high resolution (N = 5 X 10 7 ; high 
enough to obviate the need for an additional zoom-in at < lOpc scales). We 
compare a model with very weak feedback (left; q eoli =0.02 giving an effec- 
tive ISM sound speed c c ff ~ 10 km s~ ' ) to one with strong feedback (right; 
<7eos = 0.25, orc e ff~ 50km s — '). The former shows more sensitivity to res- 
olution because of sub-structure that forms in the modes; the latter gives a 
more smooth ISM so converges more quickly. At the order-of-magnitude 
level, however, the mode amplitudes and inflow rates are similar. 



gas-dominated systems; however, at least up to gas fractions ~ 0.5 
seen here, the inflow rates remain well-described by our formula- 
tion. And the Figures include both simulations with initial condi- 
tions set to have "smooth" structure (no initial modes) and those 
with a wide spectrum of initial modes inherited from the structure 
of their "parent" simulations. More than Figure [T] Figures [2] and 
HI demonstrate that our analytic model even reproduces the sign of 
M (inflow vs. outflow) reasonably well. Note also that the m — 1 
mode becomes increasingly important relative to the m = 2 mode 
at small radii (Fig. [3}. This is not only clear visually in the images, 
but also quantitatively in the mode amplitudes rl At the smaller radii 
shown in Figures[2]and[3] the modes are more often in the marginal 
orbit-crossing limit, with |£| ~ 0.1 — a few, rather than C, 3> 1; this 
highlights the importance of our results in § |2.3.3| which show that 
shocks and rapid inflow can occur even where \C\ < 1 (although not 

forC<l). 

Figure H] briefly illustrates how the simulation properties de- 
pend on the resolution and uncertain microphysics of the numer- 
ical models. We show the mode amplitudes and |£(/J)| (f° r sim- 
plicity restricting to just m — 2), and inflow rates as a function 



' The reasons for this are discussed at length in |Hopkins & Quataerq 
(2010a) and Hopkins & Quataert (2010b), but in general owe to the fact 
that k rj Q in a quasi-Keplerian potential; "slow" m = 1 modes are therefore 
generic (though not unique) to small scales near a BH (see Tremaine 2001 
for stellar-dominated systems, orlOstriker et al. 11992 for gaseous disks). 



of radius, for simulations of the central -~kpc (the range shown 
in Fig. Em. We consider otherwise identical initial conditions but 
with different resolution. We also consider two different extremes 
in the parameterization of the ISM gas equation of state, encap- 
sulated in the parameter g cos (see Hopkins & Quataert 2010a for 
details). A low value of g CO s corresponds to an ISM with efficient 
cooling and without much turbulent or thermal pressure support 
(here g C os = 0.02 corresponds to an effective sub-grid sound speed 
of just ~ lOkms"'), and a high value corresponds to an ISM with 
strong pressure support (here q eos — 0.25 corresponds to effective 
sound-speeds ~ 50kms -1 ). For each, we compare simulations at 
three resolution levels: a low-resolution run, for which the mass 
and force resolution inside of ~ 1 kpc are comparable to what can 
be achieved directly (without any "zoom-in" or "re-simulation") 
in the parent, galaxy-scale simulations; a standard-resolution run 
(the characteristic resolution of our re-simulations on this scale); 
and a very high-resolution run (the "ultra-high" resolution runs in 
Hopkins & Quataert 2010a) with ~ 0.3 pc force resolution, small 
enough to resolve the same scales as our nuclear-scale simulations 
without the need for an additional re-simulation iteration. Even at 
low resolution, the major qualitative features are in place and the 
inflow rates are order-of-magnitude accurate. But in the case with 
weak assumed feedback, the ISM is subject to much greater re- 
solved fragmentation and clumpy star formation (see Appendix B 
in Hopkins & Quataert 2010a for more discussion and examples). 
This leads to inflow rates and mode structures that are more highly 
variable in position and time, and the exact details (depending on 
how clumps collapse locally and turn into stars) are more sensitive 
to resolution. The case with stronger implicit feedback yields an 
ISM with less explicit substructure, which in turn allows rapid con- 
vergence with resolution. But in both low and high-feedback cases, 
the average behavior can still be reasonably approximated with the 
analytic prescriptions (accounting for the implicit, sub-grid ISM 
structure assumed in the numerical models). Because of the uncer- 
tain details of star formation and ISM physics in the simulations, 
we do not think any of these runs should be taken literally as an ex- 
act description of what really happens in AGN; however, they pro- 
vide a means to test our analytic prescriptions in fully non-linear, 
chaotic systems. 

Note that on small scales, we have chosen times after the 
two BHs in a galaxy merger are assumed to merge. It is of course 
not clear if the final merger proceeds quickly or slowly on sub-pc 
scales. We do so for simplicity and generality (since isolated disks 
will not have a second black hole). But during the actual binary 
phase of a near equal-mass BH-BH merger, the mode structure and 
gravitational torques on small scales will be highly non-linear and 
feature a more complex geometry than what we show here. We can, 
however, still attempt to consider this in the context of our analytic 
framework, but with the secondary BH as the "collisionless" driver 
of non-circular gas motions and shocks, as opposed to modes in the 
stellar disk. For more discussion on the inspiral phase in these sim- 
ulations and its implications for the presence of the m = 1 modes 
discussed above, we refer to Hopkins & Quataert (2010a). 

4.3 Inflow Rates 

We now present a statistical comparison between the inflow rates in 
our numerical simulations and the predictions of inflow produced 
by stellar torques (ipl. In doing so, we shall show a number of Fig- 
ures that contain the results of many of our simulations. In these 
Figures, the critical point to focus on is less the results of any given 
simulation (which can be difficult to identify in the Figure), but 
rather the ensemble properties of all of the simulations (e.g., me- 
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Figure 5. Instantaneous gas inflow/outflow rates \M\ in major galaxy 
merger simulations versus several predictors of the inflow rate. The sim- 
ulations span a range of masses, gas fractions, sizes, and orbital parame- 
ters. Some of these parameters are explicitly labeled (e.g., color denotes gas 
fraction, open (filled) points denote initial disks with structural parameters 
similar to local z ~ (high-redshift) galaxies). We calculate \M\ through a 
radius ss 5 e (where e is the simulation resolution limit; results are similar 
for radii R ~ 0. 1 — lkpc), at random times uniformly sampling between 
first passage and 0.5 Gyr after peak activity. Top: Simulation results com- 
pared to the dimensional scaling M gas (< R) Cl(R). The scatter is very large, 
~ 1.5 — 2dex. Middle: Simulation results compared to the simple estimate 
of inflow induced by strong orbit crossing in the gas (eq. |19( . Predicted and 
simulated inflow rates agree well, with a scatter of ~ 0.5 dex. Bottom: Sim- 
ulation results compared to the full scaling derived for the strong-torque 
regime (eq.|26). The scatter is reduced to < 0.3 dex. 



dian, scatter, or systematic trend relative to analytic predictions). In 
what follows, we take F(f ) = 1 for all but the m = 1 modes at < 10 
pc, for which we take F(£) = 1/2; these choices are motivated by 
the typical values of C, shown in Figures fTp] 

Figure BJ compares the instantaneous inflow rate in major 
galaxy merger simulations to several analytic predictors of the in- 



flow rate: the simplest possible dimensional expectation, M gas (< 
R)Q(R) (top panel); the order of magnitude prediction of inflow 
produced by stellar torques (middle panel; equation |19); and the 
more detailed prediction of our analytic model (bottom panel; equa- 
tion [36] from j |2.3.4} . In all three panels, we measure the inflow 
rate through a radius of 5 e (where e is the gravitational softening 
length of the simulation, roughly chosen so that 5 e ~ _R c ff/10), al- 
though we find similar results at other radii. We calculate M for 
each simulation at ~ 10 8 yr intervals, from just before first pas- 
sage until « 0.5 Gyr after the peak of activity (when the galaxies 
fully coalesce), at which point significant inflow has ceased. The 
instantaneous inflow rate is averaged over each time interval be- 
tween points to average over both the chaotic motions of the gas 
and noise introduced by the finite simulation resolution. The mode 
amplitude and other properties of the gas needed in equations [19] 
and[36]are calculated from the simulations themselves (at the same 
radius as the measured inflow rate) as in j ]4,2| 

Figure BJ (top panel) shows that although there is a correla- 
tion between the true inflow rate and M gas (< R) Q(R), the scatter is 
very large ~ 1 — 2 dex. Note that the large vertical scatter here (par- 
ticularly that above the prediction) does not imply material falling 
in on much less than the dynamical time; instead, it stems from a 
combination of time variability and scatter in the x-axis. The bot- 
tom panels of Figure BJ show that the analytic predictions in the 
strong torque regime, equations [19] and [36] work remarkably well, 
predicting the instantaneous inflow rates with a scatter of less than 
~ 2— 3. Moreover, there are no systematic trends with simulation 
properties such as the gas fraction, galaxy size, orbital parameters, 
or ISM equation of state model. 

Figure [6] compares the inflow rate to equation [19] (the order 
of magnitude strong torque M) for our other sets of simulations: 
minor mergers, isolated bar-unstable disks, intermediate (sub-kpc) 
scale re-simulations of galactic nuclei, and nuclear-scale (< lOpc) 
re-simulations of disks around BHs. In each case we again mea- 
sure the inflow rate at ~ 5e in the simulations and smooth the re- 
sults in time. The strong torque analytic model works well in all 
cases, despite a very large dynamic range in mass, spatial scale, 
and inflow rate, and despite the fact that different, independent 
non-axisymmetric modes control angular momentum loss at dif- 
ferent radii. In addition, some of the simulations (particularly the 
mergers) are clearly in the strong orbit-crossing regime (£ ^> 1), 
while others (particularly the < lOpc circum-BH disks) are in the 
marginal orbit-crossing regime (C, < 1) in which there are fewer (if 
any) formal orbit crossings. The same inflow prediction nonethe- 
less works well in both limits. Gas fractions are explicitly labeled 
by color; although we expect our assumptions to break down some- 
where around / gas ~ 1, it is immediately clear from Figure BJ that 
up to /g as ~ 0.5, there is no systematic deviation from our scalingr] 
Other parameters varied from simulation to simulation are gener- 
ally not shown, but we have checked to see that there is no clear 
systematic difference between our predictions and the results of a 
given subset of the simulations. For example, the effective equa- 
tion of state of the gas is varied widely at each scale as discussed 
in § |4. 1| but there is no statistically significant offset of the inflow 
rates with respect to the predictions in each case (with the impor- 
tant caveat that we have no simulations with a true non-sub grid 



2 We emphasize, however, that at gas fractions this high, the sub-grid 
ISM model of the simulations becomes particularly suspect (as small-scale 
clumping/fragmentation is not explicitly treated). It is safer to say that our 
results are robust to variations in / gas for / gas <g 1 . 
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Figure 6. Simulated vs. predicted inflow rates, as in Figure[5] for minor mergers, isolated galactic disks, and re-simulations of galactic nuclei; the inflow rates 
are measured at a radius of 5 e (as Figure[5j in each case, with the values of this radius in parentheses in each panel. Colors denote the gas fraction, as in Fig. [5] 
while symbol types denote some of the key parameters varied in the ensemble of simulations. In all panels, the predicted inflow rate (x-axis) is the simple 
strong orbit crossing estimate in equation [T9] Top Left: Galaxy-galaxy minor mergers. Top Right: Isolated, bar-unstable galactic disks. The parameters varied in 
the simulations include the initial disk-to-total mass ratio D/T at the effective radius and the stellar feedback prescription (g C os, where < 0.25 denotes modest 
feedback and > 0.5 denotes strong feedback; see Hopkins & Quataert 2010ai. Bottom Left: Intermediate-scale re-simulations of both mergers and isolated 
unstable systems (e.g., Fig. [2j. The key parameter regulating the formation of secondary bars is the disk-to-total mass ratio. Bottom Right: Nuclear-scale 
re-simulations of the intermediate-scale systems (e.g., Fig. 13} , where the lopsided m = 1 disk forms and regulates the gas inflow. 



model). Likewise, both in isolated disks and nuclear-scale simula- 
tions, the model of star formation has been varied in a large subset 
of simulations; the resulting points are indistinguishable. Making 
star formation less homogeneous by raising the star formation den- 
sity threshold does not appear to change our conclusions, but this 
is again with the caveat that the simulations do no explicitly model 
individual GMCs. Changes to the ISM equation of state or star for- 
mation law do affect various global properties of the results - the 
mass distribution, gas fractions, etc - but these are explicitly in- 
cluded in our formulation of M\ the important point is that they 
introduce no significant deviations from our predicted scalings. Pa- 
rameters such as the initial simulation redshift enter only indirectly, 
changing disk structural parameters, and Figure B] shows explicitly 
that they make little difference to our conclusions. 

One apparent outlier in Figure [6] is the "spray" of several 
(blue/black) points at large true M, but small predicted M, in the 
nuclear-scale simulations (bottom right). This is a system in which 
the gas disk undergoes a large-scale fragmentation and a gas-rich 
clump happens to fall on a nearly radial orbit towards the BH. This 
is a special case. Because of the high density of points in Figure 
[6] this small number of outliers makes it appear as if the scatter in 



this panel in much larger - in fact, the correlation is just as sig- 
nificant as in each of the other panels (linear correlation coefficient 
0.87, or significance > 10cr), and the intrinsic scatter is nearly iden- 
tical (lognormal 0.37 dex, as compared to 0.33 dex for the points 
plotted in other panels). Recall as well that in many of the simula- 
tions shown here, the gas motions include a significant contribution 
from resolved 'turbulent' motions and the gas mass distribution can 
be highly inhomogeneous/clumpy; we nonetheless still find good 
agreement with the predictions of the strong torque analytic model. 
Figure U\ shows the same simulation results as Figure [6] but 
compares them to the predicted inflow in the weak-torque regime, 
equations Il3| & Il7| These scalings do fine at high M, but fail 
to reproduce the numerical results at lower M and are clearly 
systematically discrepant. If we use the results of |Lynden-Bell| 
& Kalnajs (1972) for direct angular momentum transport by the 
non-axisymmetric modes in the gas (also second-order in \a\; see 
eq. [T8} , we obtain a similar, but more severe, discrepancy. This 
is because the inflow rate induced by the direct transport by non- 
axisymmetric modes (eq. [T8} is smaller than the weak-torque in- 



flow produced by the stars on the gas (eqs U3\ & 
factor of ~ Ega S /Etot (see {2.2 1. Overall, Figure 



17} by at least a 
7Thighlights that 
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Figure 7. Simulated vs. predicted inflow rates, as in Figure [6] but for the 
linear weak-torque, non-orbit crossing predictions (eq. |17| with a(A) ~ 1). 
This prediction works reasonably well at the highest M, because the second- 
order formula in \a\ approaches that for strong torques (in Fig. [6l when 
\a\ — > 1. However, over the majority of the dynamic range in M the weak- 
torque prediction significantly under-predicts the true inflow rate. Note that 
the axes here are the same as in Figure [6] As a result many of the low 
simulated M points (y-axis) do not show up in this Figure because they 
have very small predicted M (x-axis) and are thus off the panel to the left. 



for the entire mass and radial range we have simulated, the domi- 
nant torque is well-modeled as originating from shocks induced in 
the gas by the stars. 

5 BRIDGING THE GAP: FROM GALAXY TO BH 

The literature contains a number of estimates of the accretion rate 
onto a central massive BH and/or through a given radial annulus, 
none of which properly account for the physics we have elucidated 
here. Here we briefly review several of these estimates (j |5.1[ l and 
then present a new analytic model for the inflow rate onto a central 
BH given the conditions at larger radii in its host galaxy (j |5.2} . 

5.1 Previous Accretion Models 

In semi-analytic and numerical models that lack the resolution to 
study gas transport from ~ 0. 1 — 1 kpc to < 0. 1 pc, simple subgrid 
prescriptions are necessarily used to estimate the accretion rate onto 
a central BH. One common estimate of the BH accretion rate is the 
Bondi accretion rate for spherical accretion (e.g., |Di Matteo et al.| 
|2005l >: 



M B 



= 47TQB 



g 2 mI 



BH P» 



(44) 



where c s and p„ as are the ambient sound speed and gas density, re- 
spectively, and lvb is a constant that depends on the mass profile and 
gas equation of state; as is sometimes taken to be as large as ~ 100 
(Springel et al. 2005). Although perhaps appropriate for pressure- 
supported gas (e.g., in clusters of galaxies), the Bondi rate is not 
well-motivated for the cold, nearly Keplerian, gas that dominates 
the mass in disk galaxies. Indeed, we shall see that equation [44] is 
not a good approximation to our numerical results. 

For rotationally supported gas, the inflow rate is determined 
by the rate of angular momentum transport. If angular momentum 
transport in galaxies were primarily produced by local stresses in 
a thin disk, the effective kinematic viscosity would be ~ ay c s h, 
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Figure 8. Comparison of simulated and predicted gas inflow rates as a func- 
tion of radius. The simulations are shown during the strong inflow phases, 
as described in Figure [3] Top: The inflow rate in the simulations at each 
radius (dM/df) compared to the simple estimate for strong gravitational 
torques (eq. |19[ . Each solid line denotes a different simulation, with dot- 
ted lines showing the radii near the boundaries of our re-simulations, where 
the exact profile shape is suspect; colors correspond to initial gas fractions 
as in Figure [5] Middle: Simulated inflow rate versus the scaling expected 
for accretion via local viscous stresses (eq. |45). Bottom: Simulated inflow 
rate versus a Bondi-Hoyle-like estimate (eq. |44| with Mbh — ► Af enc (< R)). 
The Bondi-Hoyle and viscous accretion rate formulations do not accurately 
capture the dominant physics in the numerical simulations, while the strong 
gravitational torque estimate does quite well. 



yielding an inflow rate 

M v i SC ous = 37T av c s E gas h 



(45) 



where h ~ c s /Q is the disk's vertical thickness and ay is the di- 
mensionless viscosity parameter (assumed spatially constant). 

Both equations |44| and[45] depend quite differently on the am- 
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Figure 9. Comparison of the BH inflow rate in our simulations with the 
viscous (left; ay = 0.1) and Bondi (right) accretion rate estimates eval- 
uated at two large radii: lOOpc (top) and lkpc (bottom). Solid diagonal 
line shows equality, dotted lines ± 1 dex. Colors denote gas fraction (inside 
each radius) as in FigureB] The viscous and Bondi accretion rate estimates 
both fare poorly, with very large scatter and a large systematic offset for 
the Bondi estimate. In the simulations, the BH "accretion rate" is estimated 
using the inflow rate at < 0. 1 pc in nuclear-scale re-simulations that have 
initial conditions taken from the results of intermediate-scale or galactic- 
scale simulations with the appropriate conditions at ~ lOOpc or ~ 1 kpc, 
respectively. 



bient gas properties than the inflow rate produced by stellar gravi- 
tational torques and orbit crossing (eq.|19|&|36^. To highlight this, 
Figure [8] shows the ratio of the true inflow rate through a given an- 
nulus R in our simulations to the viscous accretion rate (eq,|45|l and 
the order of magnitude gravitational torque inflow rate from equa- 
tion|19| Because we show the results as a function of scale, all simu- 
lations on various scales are shown together (but each is only shown 
at the well-resolved and self-consistently modeled radii appropriate 
to it). We also compare the simulation results to a generalization 
of the Bondi-Hoyle formula appropriate for larger radii in a galaxy, 
with Mbh — > M cnc , the total enclosed mass within a radius. This is 
strictly valid for inflow through an annulus R only if M en c changes 
relatively slowly with R so that the spherical accretion sonic radius 
is at > R. We include it to emphasize how poorly the simple Bondi- 
Hoyle scaling does at describing the radial variation of M in galac- 
tic simulations (for discussion of how Bondi rates change when ap- 
plied outside the Bondi radius, see Colpi et al. 2007{ |Volonteri et aT| 
2008). The results in FigurelHlare for randomly chosen times near 
the peak of activity in each of the simulations from Figures [5|6] (the 
re-simulations from Hopkins & Quataert 2010al. Note in particular 
that Figure[8]shows the results of three different sets of simulations, 
with the intermediate and nuclear-scale simulations using initial 
conditions drawn from the results of the larger-scale simulations; 
the boundaries between the different sets of simulations are shown 
with dotted lines. Runs shown also include systematic numerical 
resolution studies discussed in Hopkins & Quataert (2010a). 

Figure[8]shows that the gravitational torque model works well 
at all radii, with < 0.3 dex scatter and without significant system- 
atic deviations over the range of simulations we have carried out. 
By contrast, the viscous and Bondi accretion rate fare poorly by 
comparison. Not only is there factor of ~ 10 scatter, but there are 
very large systematics as a function of radius (and in some cases, 
as a function of model parameters such as gas fraction and/or ISM 
equation of state parameterization g e os). 



Figure [9] demonstrates a similar point by comparing differ- 
ent predictors of the BH accretion rate given the conditions in 
the galaxy at large radii. Specifically, we compare the viscous or 
Bondi-Hoyle accretion rate evaluated at a large radius ^o (as would 
be done using such prescriptions as sub-grid models) to the true in- 
flow rate near the accretion disk. The latter is taken here to be the 
inflow rate at < 0.1 pc in our re-simulations that use initial condi- 
tions drawn from the galactic-scale conditions at Ro (as in Fig.|8lr| 
For the viscous accretion rate, we use av =0.1, although we con- 
sider this normalization to be arbitrary. For the Bondi-Hoyle accre- 
tion rate, we use qb = 1 , which assumes that the density and sound 
speed are relatively constant inside Ro. 

Figure[9]shows that neither the Bondi nor viscous models pre- 
dict the actual inflow rates at small radii. For the Bondi rate, there 
is essentially no correlation between the true and predicted inflow 
rates. For the viscous rate, otv ~ 0. 1 gets the median M approxi- 
mately correct, but there is factor > 10 scatter and only a weak cor- 
relation between the true and predicted inflow ratesnThe predicted 
inflow rate also depends systematically on the radius at which it 
is evaluated (as in Fig. [8](. Variations in the sound speed c s drive 
much of the horizontal scatter in the predicted Mbh in Figure M at 
a given true Mbh- This is because the Bondi and viscous models 
both depend strongly on c s (eqs.|44|&|45), whereas the true inflow 
rate is only a weak function of sound speed (at least for the stellar- 
dominated systems that are the particular focus here). 

5.2 A Gravitational Torque Model 

In this section, we use the analytic results from § [2] to derive a 
physically motivated estimate of the inflow rate M at a given ra- 
dius in a galaxy, which approximates the numerical simulations de- 
scribed in § [4] Our goal is to develop a model of the inflow rate 
that can be extrapolated over a large range of spatial scales, both to 
provide physical insight and to approximate the BH accretion rate 
on small scales from properties estimated on galactic or ~ lOOpc 
scales (e.g., in cosmological simulations or semi-analytic models). 

We assume that the galaxy properties are known at a radius 
Ro, outside the BH radius of influence. We also assume that the 
stellar disk and bulge have approximately power-law surface den- 
sity profiles, i.e., £<*,* oc R -71 ' and £(, oc R~ Vb . This is reasonably 
accurate in our simulations, so long as Ro < R e (the galaxy effec- 
tive radius), and it allows us to make analytic progress. We will 
typically consider r\b ~ rjd as well, i.e., that the bulge and disk have 
similar density profiles. 

As argued in ip] the gas in galaxies reaches a quasi steady- 
state in which its radial distribution is set by a competition between 
inflow to smaller radii and star formation (Thompson et al. 2005 1: 

^-M M0 w~27r#E, . (46) 

oR 



3 In a simple alpha-disk prescription, the inflow rate onto the outer edge 
of the disk corresponds to the actual accretion rate onto the BH. In prac- 
tice, there may be a complex transition to this inner disk, for which some 
intermediate-scale model (e.g. Power et al. 2010) should be applied to map 
the efficiency of inflow through to the viscous disk. 

4 It may appear surprising that the local viscous estimate is even accurate to 
a factor of ~ 1 00. This is because the sound speed in our simulations is 3> 
the true thermal speed of neutral molecular or atomic gas (which dominates 
the mass); using the latter for c, would lead to a much smaller M v iscous- Our 
viscous estimate is thus more akin to an estimate of the inflow rate produced 
by turbulence generated by stellar feedback or strictly local, sub-sonic grav- 
itational instabilities jGammie|2001| jLodato & R ice 2004 Durisen et al.| 
[2007) . 
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The qualitative behavior of solutions to equation Bm depends on 
the ratio of the star formation time U to the inflow time tmf. For 
ftnf ^ U, the solution is Minflow — const ~^> M,, while if Anf S> r», 
Mi n fl ow declines precipitously at smaller radii as gas is consumed 
by star formation. The observed star formation times in galaxies 
are ~ 10 — 100 dynamical times, which is comparable to the inflow 
time for a ~ 0.01 — 0.1. This is similar to the mode amplitudes 
in our simulations and thus we expect t- m ( ~ f* . A reasonable ana- 
lytic approximation to the numerical solutions of equation «[46jl in 
this limit is given by taking Mi n fl ow ~ M* ~ nR'T,, . More precisely, 
this approximate solution is valid when it itself implies Mi n fl ow oc R p 
with p > 0, so that Mi„fl ow declines at small radii (as it must in the 
correct solution). This is not guaranteed for all r\d,r)b,r\K, etc., but is 
satisfied for the parameters we utilize below. The exact normaliza- 
tion will depend on details of the star formation law beyond just the 
power-law scaling (e.g. small-scale structure not treated explicitly 
in the simulations), so we consider it somewhat uncertain. 

In the potential of the galaxy, the dominant mode of angular 
momentum transport is an m = 2 bar mode. To estimate Mj n n ow , we 



use equation 36 with the WKB approximation for |$i 



l*i 



,2nGT, d R\kR\ 



(47) 



We will eventually take kR ~ 1 because the most important modes 
are global. Since we are focusing on small scales in the galaxy 
(< 0. 1 — 1 kpc), we would like to know the saturation mode am- 
plitude |a| max appropriate for self-excited modes (rather than the 
mode amplitude induced directly by an external perturber). This is 
in general a complex function of non-linear processes such as heat- 
ing by the mode, fragmentation, inflow, and momentum exchange 
with the halo. Hopkins & Quataert (2010a) show the saturation am- 
plitude |a| m ax in simulations on this scale as a function of various 
properties of the galaxy (their Figure 12); the results can be reason- 
ably approximated by |a|max ~ aofd where f d = 7rGE<j/f2 R (a 
M d {R)/M cnQ {R) and ao ~ 0.3. Given this, equation [46] can now be 
approximated as 






(48) 



mao\2w +Vs s -3m\(2+vs t ) z F(() 

a d = ', 7^T. rTTTTl (49) 



4n(2 + v n )\kR\ 



v u \ 



<91m< 
a In/? 



(50) 



where Ea: and tg are the parameters of the Schmidt-Kennicutt law 
(eq.[39|, E, = E rf + E b , and we take n 2 (R) « GM enc (< R)R~ 3 . Of 
course, the dominant term need not be just a bar, but the important 
point is that the dimensional scaling is the same for each - in the 
linear regime, modes are independent so each mode enters as ~ 
ma,„F{C,m) and the relevant normalization reflects the sum over the 
spectrum of modes (in simulations, this typically does not differ 
from the single-mode normalization by more than a factor of a few). 
Near where BH dominates the potential, an m = 1 mode 
often dominates the transport. Again, in real systems there can 
be a rich mode structure, but we focus on this case for conve- 
nience and tractability. For a power-law mass profile, M(< R) = 
M mc (< Ro) (R/Ro) v , it is straightforward to solve for the ra- 
dius Rbh at which the BH dominates the potential, i.e., where 
M(< Rbh) = Mbh- The inflow rate in towards the BH from this ra- 
dius then follows from equation|36] but now using the coefficients 
evaluated for an m — 1 mode in a roughly Keplerian potential. Be- 
cause the modes here are global, we assume that $i(< _/?bh) = 
a GMd {Rbh ) /Rbh in the BH potential (the local/WKB term is small 



as R — > 0). The approximate analytic solution of equation l |46[ l is 
again M lnfl0 w(^) ^ M*(R), which implies T, g (R < Rbh) oc R' vbh , 
with r/BH = 1/2(t)k — 1). Specifically, in the limit of large gas sup- 
ply, this leads to 



E* 



VK-l 



' Qbh 



U- 



^te(^) /<(**=) (5D 



tf BH '* V^b H ; 



«BH =fllF((l 



(52) 



3 r\ K — 4 2 + vy. k 
2ty( Vk -1) 2~ 

where a\ ~ 0.2 is defined by |aj max = |$i(#bh)|/$ ~ a\ fd{Rmi) 
and the numerical value is based on the simulations of modes in the 
BH potential from Hopkins & Quataert (2010a). Together with the 
solution for Rbh using equation[48] this in turn gives the following 
BH inflow rate at a radius R: 



M = aZ K Rl£l(Ro 
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BH 



(53) 

(54) 
(55) 
(56) 
(57) 



Equation[53]is an estimate of the BH inflow rate in the limit of 
large gas supply (because we have used the Mi n n ow ~ M, solution 
to eq.|46[). If the gas supply at Rbh is less than a critical amount the 
inflow rate will be smaller than given by equation [53] The critical 
gas mass M gas ,c which defines the large gas supply limit is found by 
integrating E gas in equation |5T] from R — to R = ^?bh- Recall as 
well that the intermediate-scale physics described by equation [48] 
itself sets the amount of gas that reaches Rbh- The gas mass sup- 
plied to this radius, M gas ,bh, can be estimated by integrating E gas as 



implied by equation 48 over the same range. For M gas .bh > M gasc , 
equation[53]for the BH inflow rate is reasonably accurate while for 
Mgas.bh < Afgas, c it overestimates the inflow rate. To quantify the 



46 



suppression in M for M gas .bh < M gas . c , one can solve equation 
as a boundary value problem with a given gas mass M gas .bh at Rbh 
(as, e.g., [Thompson et al. 2005 did using somewhat different as- 
sumptions about the inflow and star formation physics). For low 
gas masses the solution is an inflow rate that is approximately in- 
dependent of radius and oc M gas .bh, implying that the suppression of 
M is given by 



M^M 



M gas . c 

-M gas .bh 



M, 



gas,c 



M, 



= C*ri 



gas.bh 



4riK — 5 — T)d faBH\' 1 K- 1 



4r) K -5 



<•■>■/ 



(58) 
(59) 
(60) 



Following the same logic, a similar suppression factor in M 
should appear if the gas mass available at Ro is too small. Specif- 
ically, the solution for the gas surface density profile between Rbh 
and Ro in equation [48] is valid only in the limit of large mass sup- 
ply. It requires that the available gas mass at Ro, M gas (/?o), exceed 
a critical value M gaSj i that can be found by by integrating E gas in 
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equationB8]from R — to R = Ro. The resulting suppression factor 
is then 



M^tM 



1 + 



M eas 



M gas (tf ) 



M a . 



4tt (t7at — 1 ) 



(n^o)^)"*- 



(61) 

(62) 
(63) 



5 4^-^-5 " 

We now combine terms and evaluate the numerical coeffi- 
cients and exponents in equations [53] and [58] We adopt star for- 
mation parameters consistent with those measured by Bouche et al. 



(2007): Y,k = lO 8 M kpc _2 , t K = 0.41 x lO'yr" 1 , and t)k = 7/4 



We find reasonable consistency with our simulation results and the 
profiles of observed cusp ellipticals by taking rji, « T)d w 1/2 over 
this radial range (see e.g. Lauer et al. 2007. Hopkins et al. 2009a). 
We assume that the R > Rbh mode is m = 2, while the R < Rbh 
mode is m = 1, and that both are global, so that \kR\,C, ~ 1. We 
emphasize that all of these are gross simplifications and shouldn't 
be taken to be universally true - however they provide a starting 
point for the order of magnitude normalization of our derived in- 
flow rates. 

Finally, equation {53} should be evaluated at the radius W acc 
where accretion itself maintains stability to local gravitational per- 
turbations, i.e., g> 1. Inside this radius, star formation ceases, 
and M ~ constant. The value of Race in fact depends on the mecha- 
nism of angular momentum transport. For local angular momentum 
transport, i? acc ~ 0.01 pc, while for larger-scale torques, R xc ~ 0.1 
pc ([Goodman 2003). Our calculations favor the latter physics, so 
we take R dcc ~ O.OI^bh, which corresponds to i? acc ~ 0.1 pc for 
typical parameters (it makes little difference, to within the scatter 
and accuracy of our formula, if we take R ia: to be a fixed physical 
radius instead). 

With these assumptions we find that our analytic gravitational 
torque inflow model predicts a BH inflow rate of 
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where a{rjK) is a dimensionless constant discussed below. In a sim- 
pler form, equation[64]becomes 



MBH.grav 

MQyr" 1 
with 



3/2, 



<*(m) fy-M^M^R-^ (1 +ft/fm) (65) 



/o~0.31/j 



M d (R ) 
10" M 
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-1/3 



f m (Ro)=M m (Ro)/M d (R ) 



(66) 
(67) 



Recall that fy in equations [64] and [65] is the total (stellar and gas) 
disk mass fraction evaluated at Ro- Note also that these predictions 
for MBH.grav can, under a variety of circumstances, exceed the Ed- 
dington accretion rate. How much above Eddington the BH can in 
fact accrete depends on the uncertain physics of radiation-pressure 
dominated, super-Eddington accretion; much of the mass may in 
fact be unbound at radii smaller than we consider here. This is less 
likely to be a concern for MBH.grav < Afeid- 
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Figure 10. Comparison of the BH inflow rate (into < 0. 1 pc) in our simula- 
tions with the predictions of our analytic model of inflow driven by stellar 
gravitational torques (eq. |65| with a(rnc) = 5). In each panel, the proper- 
ties used to infer Mbh are measured at a different radius Ro, from 30pc 
to 1 kpc (each from simulations modeling the appropriate scales). Colors 
denote gas fraction (inside the given radius) as in Figure B] The analytic 
prediction works reasonably well and the agreement is much better than for 
the Bondi or viscous predictions at all radii (see Fig. [51. There remains, 
however, order-of-magnitude scatter - mostly owing to time variability in 
the BH inflow rate that should be included in models. 



Despite the complex expressions leading to the final result for 
MBH.grav, the power-law exponents in equations 64 and 



65 



are rea- 
sonably insensitive to the values of r/j and t\k over the plausible 
range (0 < r/j < 1, 1.4 < t\k < 1.75). The normalization, however, 
is fairly sensitive, especially to r/ K , ranging from (a 1 Mq yr to ss 
18Mq yr for t\k = 7/4 to 3/2, respectively (for r/j — 1/2). This is 
one reason why we explicitly include the constant a(r]x) ~ 1 — 10 
in equations [64] and [65] The dependence on r\ K is due to the fact 
that a larger tjk implies a higher star formation rate in dense gas 
and thus a lower BH inflow rate. Star formation on the small scales 
of interest here is uncertain - we choose to parameterize that in 
a so that the model has the advantage that it can be rescaled in a 
straightforward manner for many different models of star formation 
(via an appropriate choice of t]k or r^-). We have also made the over- 
simplified assumption of just one mode at each scale, when the real 
behavior reflects a sum over a mode spectrum; this does not change 
the dimensional scalings, but does enter the normalization a. The 
exact values we adopt here are therefore unlikely to be universally 
valid, but are convenient, observationally motivated, and match the 
choice in most of our simulations. 

Figure[l0]compares equation[65] with a(»jjt) ~ 5 and with the 
relevant quantities evaluated at different Ro in the simulations, to 
the actual inflow rates into < 0. 1 pc. There is, of course, very large 
scatter - even with full knowledge of the mass profile and mode 
amplitudes, the inflow rates are chaotic and time- variable. And un- 
surprisingly, the accuracy of equation[65]decreases for larger values 
of Ro- Moreover, because the profile shapes and mode structure do 
not exactly match what we have assumed at all radii, the appro- 
priate normalization of MBH.grav depends weakly on Ro. In spite of 
its simplicity, however, equation [65] captures the broad trends in 
the BH inflow rate, and is a much more accurate approximation 
than either the viscous accretion rate or Bondi-Hoyle accretion rate 
(compare Figures[T0]and[9]l. 

As a sanity check on the predictions of equation[65] we briefly 
use it to estimate the inflow rate for two well studied systems: 
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the Milky Way and NGC 1068. For the Milky Way, which has 
M BH ~4x 10 6 M q ( JGhez et al.||2008| and /„ , ~ 0.4, f d ~ 0.2 
and M ea5 (/?o) ~ 3 x_L0'M Q 
ris|l996i, equation 



65 



at R ■ 
predicts MBH,grav 



100 pc jSerabyn & Mar 

J M yr 



Al- 
this 



3x 10 
though much larger than the current inflow rate onto Sgr A 
is not unreasonable for a time averaged inflow rate given that BHs 
of ~ 10 7 M© are growing significantly in mass in the local uni- 
verse peckman et~aLl[2004| For NGC 1068, a nearby Typ e 2 
Seyfert with M BH ^ 1.5 x 10'Mr JGreenhill & Gwinn|l997t and 
L ~ 10 45 ergss~' ~ Lz&a 'Pier et al.||!994f, the gas mass within 



Ro ~ 100 pc is ~ 3 x 10 Ma whi le the total dynamical mass is 
~ 3 x 10 8 M© JHelfer& Blitz 1995). The total disk mass within this 
radius is uncertain but Davies et al. ( 2007 1 find evidence for a stellar 
10 8 M Q within ~ 100 pc (or less). Adopting / gas ~ 1/3 

0.03 M yr -1 ' 



disk with 

and fd ~ 0.4, equation 



65 



predicts MsH.grav ~ 0.0:5 M© yr ' compa- 
rable to the Eddington accretion rate ~ 0. 1 Mq yr - ' for NGC 1068. 
There is, of course, no reason to expect that our estimate should ac- 
curately predict the inflow rate in a given system at a given time 
(as indeed it does not for the Milky Way!). It is, however, reassur- 
ing that our estimates of M for the Milky Way and NGC 1068 are 
observationally plausible in a time-averaged sense. 

It is worth highlighting how and why our inflow rate predic- 
tion (eq.[65j depends on several key physical parameters: 

• fj\ The inflow rate decreases rapidly with decreasing disk 
fraction (i.e., increasing bulge to disk ratio). This is expected for 
inflow driven by global gravitational instability in the disk. 

• Mbh: The absolute inflow rate is not a strong function of Mbh. 
In units of the Eddington accretion rate, the inflow rate actually 
decreases with increasing BH mass (~ M^ for r)a = 1/2). This 
is because on small scales, the BH plays the same role as the bulge, 
suppressing global modes within Rbh- 

• Mj: MeH.grav is linear in Mj, the total (gas+stellar) disk mass. 

• /gas: Provided that sufficient gas is available (/ gas > /o), the in- 
flow rates are nearly independent of the gas supply, because the gas 
densities at small scales are set by an equilibrium between inflow 
and star formation, which self-adjusts to a particular inflow rate. Of 
course, at very high / gas , our assumption of stellar-dominated disks 
breaks down, and it is not clear how robust the extrapolated model 
will be. 

• M* : The star formation rate enters into the expression for 
MBH.gmv via the values of tK, £k, and tjk from the Schmidt- 
Kennicut law (see eq. [53] for the precise dependencies). A lower 
star formation efficiency (corresponding to larger tK and/or smaller 
tjk) at small radii in gas-rich galactic nuclei, as several authors have 
suggested (Thompson et al. 2005; Begelman & Shlosman 2009 
Larson 2009}, would significantly increase the final BH accretion 
rate, making it easier to grow a central massive BH. 

• c s '. In our model, the inflow rate is independent of the sound 
speed c, and other thermal properties of the gas. This is because 
the torques and inflow rate are determined by supersonic non- 
axisymmetric motions induced by global gravitational instabilities. 
Again, caution is needed when systems have very large / gas and 
fragmentation can compete with global torques. 

Evaluating equation [65] requires the input parameters Mbh, 
M gas (Ro), Md(Ro), Mt(Ro) and a choice of Ro (ideally the smallest 
resolved radii, since Fig.[T0]shows that the analytic approximations 
here are best at smaller radii < 100 pc - but they are still far su- 
perior to the Bondi and viscous models even at ~ kpc). In analytic 
and semi-analytic models, these are all clearly-defined, and imple- 
menting equation [65] is straightforward. In numerical simulations, 
it is less clear how to define the bulge and disk separately in an on- 



the-fly manner. One possibility that may be sufficient at the level 
of accuracy of our model is to define the disk mass fraction via the 
total angular momentum of the material within Ro of the BH. 

As noted above, there is considerable variability on a range 
of timescales in our simulations of AGN fueling - see Hopkins & 
Quataert (2010a! for a more extensive discussion. Very roughly, 
we find approximately equal power in fluctuations on each dynam- 
ical time resolved in the simulation. Recently, |Levine et al.| (2010) 
show this explicitly for timescales from > 10 7 to ~ 100 yr (see their 
Figures 3-5) using AMR simulations that cover a range of radii 
similar to our simulations. Given this variability, it is unlikely that 
any BH accretion rate model could predict inflow rates on small 
scales given the properties of a galaxy on large scales with much 
better than order-of-magnitude accuracy (as in Fig.|10[>. Moreover, 
this variability can have important consequences for the effects of 
AGN feedback or the ionization of the intergalactic medium, to cite 
just two examples. We thus believe that it is important to include 
variability when modeling BH growth and AGN feedbackrj 

5.3 Implications of our Gravitational Torque Model 

Having derived the new inflow rate predictor in equation [65] and 
shown that its functional form is significantly different from the 
Bondi or viscous accretion rate estimates, we briefly speculate on 
conditions under which we do or do not expect this to make a sig- 
nificant difference relative to previous results. Because the problem 
of black hole fueling with feedback is highly nonlinear, it is entirely 
possible that the results of applying equation [65] in numerical sim- 
ulations could be more subtle than we anticipate. Thus it will be 
critical to explore this numerically as well. 

We anticipate that our new fueling model may only have lim- 
ited effects in the following circumstances: 

• Eddington-Limited Accretion: In e.g. simulations of galaxy 
mergers, the AGN passes through a significant phase in which very 
large inflows from galactic scales lead to a large gas reservoir in the 
central ~ 100 pc. The implied accretion rates for a range of models 
are comparable to, or in excess, the Eddington rate (Di Matteo et al. 
|2005||Debuhr et al.|20l0] >. This suggests that the Eddington limit is 
the rate-limiting factor in many cases, not the rate of mass supply. 

• Feedback-Regulated Accretion: Once BHs reach a certain 
mass, models that include BH feedback often lead to self-regulated 
BH growth. In this limit, both simulations and analytic arguments 
suggest that the accretion rate is determined by the balance between 
feedback and inflow, independent of the precise functional form 
of the accretion rate model (Debuhr et al. 2011). Because the BH 
gains most of its mass during this phase, we suspect that the new 
accretion rate model proposed here will not significantly change 
previous inferences about the conditions required to reproduce the 
Mbh — a relation and other BH-galaxy correlations. 



5 For example, following |Levine et al.| J2010J, in semi-analytic mod- 
els or galactic-scale simulations one could approximate the variability on 
timescales smaller than some resolved timescale to ~ ^o(^o), and down to 
a minimum variability timescale dt, using a Fourier series 



. TT i <r lnl0 
M = (M) I I exp -j — ; ft cos 



=0 



Oat 



{<TLodt)</ N 



-2ic5i 



(68) 



where cro ~ 0.5 — 0.6 is the amplitude of variations (motivated by the simu- 
lations) and the /, and 5; are a Gaussian random numbers drawn with disper- 
sion of unity and uniform random number from 0— l, respectively; N > 10 
appears sufficient for a reasonable approximation of the variability in the 
simulations. 
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There are, however, other contexts in which we suspect that 
the more physical treatment of BH fueling presented here will have 
a significant effect on the results of analytic models and galactic- 
scale simulations: 

• AGN Lightcurves: Most BHs, most of the time, accrete below 
the Eddington rate. In this limit, the predicted AGN accretion rate 
is more sensitive to the details of the accretion rate model. As a 
result, our new accretion rate predictions may significantly affect 
the AGN lightcurves and the fraction of time a BH spends at a given 
accretion rate. As emphasized by Hopkins et al. (2005 ), this will in 
turn affect the interpretation of the AGN luminosity function. 

• Early BH Growth: The Bondi accretion rate is oc M BH , which 
suppresses accretion onto low mass BH seeds at high redshift. By 
contrast, the model presented here depends only weakly on BH 
mass - low BH mass in fact allows global instabilities to form more 
readily and thus should promote rapid seed growth. 

• BH Growth Versus Galaxy Morphology: Our accretion model 
includes an explicit, strong dependence on the nuclear morphol- 
ogy of the galaxy. BH growth via secular processes will be rel- 
atively more efficient in disk-dominated galaxies, but suppressed 
in bulge-dominated galaxies, potentially requiring more violent 
events (mergers) to grow BHs. Because the morphology as a func- 
tion of radius is also important, the bulge structural properties are 
also relevant for the resulting BH accretion rate: e.g., low mass but 
sufficiently compact bulges can suppress accretion on small scales, 
despite allowing inflow to the central ~kpc. This could introduce 
differences between the predicted accretion rates in different types 
of bulges (e.g. classical or pseudo-bulges). The dependence of the 
accretion rate on galaxy morphology may end up being one of the 
most important differences between the predictions of our model 
and previous results. 



6 DISCUSSION AND CONCLUSIONS 

The dominant mechanism generating significant gas inflows from 
~ lOkpc to ~ lOOpc in galaxies has been extensively studied us- 
ing numerical simulations over the past two decades: a perturbation 
to the galaxy (e.g., a merger or disk instability) generates a non- 
axisymmetric structure. If the perturbation is sufficiently strong in 
the collisionless component of the system, it drives a correspond- 
ing perturbation in the gaseous component. This, in turn, forces 
the gas to shock, which dissipates energy and relative motion, al- 
lowing for net angular momentum loss and inflow. We have re- 
cently demonstrated that qualitatively similar physics can operate 
from ~ 0.1 — 10 pc in galactic nuclei if they also have a domi- 
nant collisionless (stellar) component, although the relevant non- 
axisymmetric mode is preferentially in — 1 instead of m = 2 as on 
larger scales (Hopkins & Quataert 2010a). The critical physics of 
stellar-induced shocks in gas has thus been seen in simulations of 
galaxy-galaxy mergers, unstable galactic disks, and secondary and 



re-distribute angular momentum via resonances iLynden-Bell & 



tertiary gravitational instabilities in galactic nuclei (Noguchi 1988 
Barnes & Hernquist 1996; Barnes 1998; Berentzen et al. 2007 
Hopkins et al. 2009b, Hopkins & Quataert 2010a). In this paper, 
we have developed an analytic framework to model this physics 
and to describe how gravitational perturbations redistribute angular 
momentum, and drive shocks and gas inflow, in galaxies. We have 
also tested these models against a large suite of numerical simula- 
tions. 

Dissipation and the co-existence of collisional (gas) and 
collisionless (stars) disky components are critical for generat- 
ing large gaseous inflow rates. Dissipationless components can 



Kalnajs|1972 ; Goldrei ch & Tremaine|1 979 , Hernquist & Weinberg 
1992| |Athanassoula 2003, Weinberg & Katz||2007| > but these are 
much less efficient than gravitationally-induced shocks, with inflow 
rates smaller by a factor of > 10 — 100. And without a significant 
fraction of the mass in a collisionless component that drives the 
gas into shocks, the gas orbits would self-adjust to eliminate any 
long-lived strong shocks. 

The canonical criterion for gas orbits to shock in the pres- 
ence of a non-axisymmetric potential perturbation is that CI — 
\dR\/dRo\ > 1, where Ri is the perturbation to a circular orbit ini- 
tially at radius Ro due to the non-axisymmetric potential (eq. |20|l. 
This criterion is, however, only a sufficient condition for shocks, 
not a necessary one. In the case of shocks near spiral arms or bars in 
a self-gravitating disk, the traditional condition on \(\ is often ade- 
quate, since |£| ^> 1 for moderate-amplitude perturbations. But this 
is not always the case, particularly for disks in a quasi-Keplerian 
potential near a BH. We have demonstrated that even when |£| < 1, 
a stellar perturbation can produce shocks in a gaseous disk with 
a finite sound speed (eq. |29|30| . Physically, the key point is that 
when the non-axisymmetric potential is stellar dominated, gas "col- 
lisions" can happen not just at strict orbit crossings, but at a finite 
separation between streamlines. The key requirements are that the 
amplitude of the non-axisymmetric stellar perturbation exceed the 
critical value given in eq. [30] and that a reasonable fraction of the 
mass is in a disky collisionless component. We have not quantified 
in detail how large the mass fraction in the collisionless component 
must be, but simple force considerations suggest > 50%. 

Given the presence of shocks in the gas, we have derived the 
resulting angular momentum loss/gain and the gas inflow/outflow 
rates in two distinct limits: very strong shocks (|£| 3> 1) and 



marginal orbit crossings (|C| < 1) (eq. 36 1. This calculation includes 



a proper treatment of the phase function for angular momentum ex- 
change between the gas and stars, i.e., when the net torque produces 
gas inflow and when it produces outflow. For typical conditions in 
galaxies, we find flow towards an inner Lindlad resonance inside of 
co-rotation and towards an outer Lindblad resonance outside of co- 
rotation. However, the behavior can be more complex depending 
on the details of the system. 

One of our key results is that the gas inflow produced by 
gravitationally-induced shocks is robust and rapid, with 



U< 



\a\M m 

'civil 



(69) 



where jaj is the fractional non-axisymmetric perturbation to the po- 
tential (eq.|19|&|36}. Note several key facts about this scaling. First, 
M is linear in \a\. In contrast, when shocks do not occur, angular 
momentum transport by non-axisymmetric perturbations is at least 



second-order in jaj (5 2.2 1. Moreover, this second-order transport is 



also suppressed by powers of kR and m where k (m) is the radial 
(azimuthal) wavenumber of the non-axisymmetric mode (see, e.g., 
Kalnajs 197 Is exact solution for the logarithmic spiral; eq |18|). 
Given these differences, we find inflow rates that are a factor of 
~ 10 — 100 larger than would traditionally be predicted by linear 
theory (for a given set of gas properties). 

A second key property of equation[69]is that the inflow rate is 
independent of the sound speed of the gas c s , i.e., of the gas thermo- 
dynamics. Physically, this is because the angular momentum trans- 
port is governed by the rate at which gas orbits are forced to shock 
via the non-axisymmetric stellar potential. Because the relevant gas 
motions are by definition supersonic and induced by a collisionless 
component, the precise gas thermodynamics is not important. This 
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property of inflow produced by gravitationally-induced shocks is 
fundamentally different from angular momentum transport by lo- 
cal viscous stresses ("a") or by linear spiral waves, both of which 
depend on CjljThe lack of an explicit dependence on the gas ther- 
modynamics, which is particularly uncertain in a multi-phase in- 
terstellar medium with stellar feedback, points to the robustness of 
our predicted inflow rates in this regime. One potential subtlety is 
that our analytic work assumes that a significant fraction of the ISM 
consists of relatively diffuse gas, rather than gas bound up in, e.g., 
star clusters; this is also true to some extent in our numerical work 
(see j ]4.1 1 for a discussion of this point). This assumption is moti- 
vated by the observed low star formation efficiency in dense gas and 
the small fraction of a molecular cloud's mass that is converted into 
stars before the cloud is disrupted (e.g., Krumholz & Tan 2007). 
But further numerical work including a more detailed treatment of 
stellar feedback and a multi-phase ISM would shed considerable 
light on whether a highly inhomogeneous ISM qualitatively mod- 
ifies the physics of gas inflow in galaxies. It is well-known, for 
example, that mode development depends significantly on gas ther- 
modynamics in pure gas disks (e.g. Gammie 2001 , Nayakshin et al. 
2007; Cossins et al. 2009). We suspect that the physics described 
here will be applicable to any collisionless plus collisional disk sys- 
tem, provided that the post-shock cooling time of the gas is short 
compared to the local dynamical time and the collisionless com- 
ponent is gravitationally dominant (well-satisfied in galaxies at the 
radii of interest). 

Another way to appreciate the importance of equation [69] 
is to compare the resulting inflow time of gas in galaxies (~ 
a| '?dyn) to the star formation timescale. The latter is observed 
to be ~ 10 - 100 fdyn (e.g., |Bouche et "aT1|2007} , Since the non- 
axisymmetric mode amplitudes produced by mergers or instabili- 
ties of self-gravitating disks can readily reach \a\ ~ 0.01 — 1 (Fig. 
12 of Hopkins & Quataert 2010al, equation[69]implies that the gas 
inflow time can be comparable to or shorter than the timescale for 
the gas to turn into stars. This is the fundamental reason that we find 
BH inflow rates ~ 1 — WMq yr _l in our simulations (sufficient 
to fuel luminous quasars) without overproducing star formation in 
galactic nuclei. Similar inflow rates have also been found in inde- 
pendent simulations by |Levine et al. (2010). By contrast, a number 

Goodman 2003; Sirko & Goodman 
Nayakshin & King 2007) assumed 



of previous calculations (e.g., 
|20T)3"1 [Thompson et al.||2005| 



either a local a model or linear transport by spiral waves, neither 
of which accurately describe our analytic or numerical results. Be- 
cause of this incorrect treatment of angular momentum transport, 
these previous calculations also significantly overestimated the in- 
flow time of gas in galaxies, and thus overestimated the fraction of 
the inflowing gas that turns into stars. In our simulations, we find 
that star formation and inflow reach an approximate equilibrium, in 
which the inflow rate is comparable to the local star formation rate 
(except at sufficiently small radii, where the inflow rate becomes in- 
dependent of radius). This balance determines the quasi-steady gas 
surface density as a function of radius (§ [3l, and is in good agree- 
ment with the numerical results from ~ 10 -3 — 10 kpc (e.g., Figs 
[2]&[3]l. Of course, star formation in these models is still "sub-grid" 
and pegged to an empirical law, so avoiding runaway local collapse 
may still require processes such as stellar feedback. 

For the purposes of analytic or semi-analytic calculations and 
galactic or cosmological numerical simulations, it would be use- 



6 More precisely, transport by linear spiral waves only depends on c s indi- 
rectly, if the system maintains Q ~ 1; see eq.|18| 



ful to have an estimate of the BH inflow rate based only on galaxy 
properties at much larger radii (generally > lOOpc). Calculations 
in the literature often assume simple models such as Eddington- 
limited accretion for a fixed timescale, spherical Bondi accretion 
( Springel et al. 2005 , Di Matteo et al. 2005 ), local viscous transport 
(e.g.,|Sirko & Goodman 2003 , Debuhr et al. 



port by spiral waves (e.g., Thompson et al. 



2010), or linear trans- 
2005). None of these 



models, however, accounts for the key physics of gravitationally- 
induced shocks in the gas. We have explicitly demonstrated that 
both spherical accretion and local viscous transport are relatively 
poor approximations to the inflow seen in our numerical simula- 
tions (Figs.[8]&[9]l. 

Using our analytic models of gravitationally-induced gas in- 
flow we have developed a more accurate predictor of the BH ac- 
cretion rate given the large-scale properties of a galaxy (eqs [64] - 
[67} . This model is highly simplified, and clearly far from exact, but 
it includes well-motivated dependencies for some of the key pa- 
rameters of the problem. For example, the inflow rate is strongly 
suppressed for larger bulge-to-disk ratios, as is generally the case 
for global gravitational instabilities. In addition, the inflow rate is 
essentially independent of the gas sound speed so long as cooling 
is efficient and most of the gas resides in a rotationally-supported 
disk. By again comparing to our suite of numerical simulations, we 
have demonstrated that this new BH inflow predictor represents a 
major improvement relative to previous models (Fig.|10|>; see j j5.3| 
for a brief discussion of some of the implications of this result. 
However, even this new BH predictor has a factor of ~ 10 scat- 
ter when compared to the simulations. This does not appear to be 
due to any systematic dependence on galaxy structural properties, 
gas fraction, etc. that is not already captured in our scalings. Rather, 
the scatter is a consequence of physical variability in the inflow rate 
for a given set of conditions in the galaxy at large radii. Indeed, be- 
cause gravitational instabilities operate on all scales, and because 
the systems of interest are chaotic, we find that there is variability 
on essentially all timescales (see also Levine et al. 2010). We be- 
lieve that it is important to include this variability in future models 
of BH growth and AGN feedback. 

A number of important regimes of gas inflow require signif- 
icant further study. We have assumed throughout this paper that 
the non-axisymmetric mode amplitudes are large enough to pro- 
duce shocks at nearly all radii. This is indeed borne out by our (and 
others') simulations of galaxy-galaxy mergers and gas-rich galac- 
tic nuclei, but it is by no means clear that this will always be the 
case. For low mode amplitudes, we do not expect shocks and so the 
angular momentum transport mechanisms studied in this paper will 
not be present. Understanding which processes dominate inflow at 
lower inflow rates will be of considerable interest for models of 
low-luminosity AGN and weak enhancements in nuclear star for- 
mation. 

The regime in which galaxies are gas-dominated is also of 
considerable interest. In this limit, it is unlikely that strong or- 
bit crossings and/or coherent long-lived shocks will occur, if the 
mode remains linear. Generically, we expect this to decrease the 
efficiency of gas inflow. This has in fact been seen and studied in 
detail in the context of galaxy-galaxy mergers (e.g., |Hopkins et aT] 
2009d). Of course, if strongly non-linear modes can be sustained, 
then our assumptions are probably not applicable and there can in 
principle be a wide range of inflow rates - a pure radial inflow at 
the free-fall velocity is a valid non-linear configuration! And in very 
gas-rich systems, local fragmentation may be much more efficient, 
generating an inhomogeneous ISM that may provide mechanisms 
for angular momentum transport that we have not considered in 
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this paper. Understanding this physics on smaller scales in galactic 
nuclei is particularly important for high-redshift galaxy formation 
and the formation of the first BHs, where it has been posited that 
they might form in direct collapse of primordial gas clouds (Begel- 
|man & Shlosman 2009 and references therein). In our simulations 
of gas dominated systems to date (Hopkins & Quataert 2010a), we 
find that the gas inflow from ~ 10 pc to ~ 0. 1 pc is significantly 
more efficient after some fraction of the gas has turned into stars. 
This is qualitatively consistent with the results from galaxy merger 
simulations, but needs to be studied in more detail (for comparison, 
see e.g. Agertz et al. 2009, Teyssier et al. 2010). 
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